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ABSTRACT 

\, 

'''^Numerical  solutions  of  supersonic  viscous  flows  are  studied  by  applying  an 
implicit  time-dependent  scheme  to  the  thin-layer  Navier-Stokes(TLNS)  equations. 
The  alternating  direction  implicit(AOI)  scheme  is  first  formulated  to  solve  transonic 
viscous  axisymmetric  flows  in  two  dimensions.  The  results  indicate  that  the  ADI 

scheme  is  not  efficient  enough  for  supersonic  viscous  calculations 

NT  "  ’  ’  -  -  -  ...  J 

Accordingly,  a  spatial  discretization  scheme  using  upwind  flux-vector  split  dif¬ 
ferencing  in  the  streamwise  direction  and  central  differencing  in  the  cross-stream 
direction  is  chosen.  Three  approximate  factorization  schemes  and  one  fully  implicit 
direct  solver  are  considered.  Of  them,  the  diagonally  dominant  ADI(DDADI)  and 
the  parabolized  ADI  are  found  to  be  much  faster  than  the  standard  ADI  procedure. 
The  optimum  CFL  number  for  the  DDADI  method  is  about  5000  and  it  provides 
competitive  convergence  with  direct  solvers.  In  terms  of  CPU  time  requirements, 
the  parabolized  ADI  procedure  is  as  fast  as  the  DDADT  method. 

st*hese  numerical  algorithms  are  applied-to  solve  supersonic  flows  through  coni¬ 
cal  and  high  expansion  ratio  contoured  nozzles  for  different  Reynolds  numbers,  wall 
temperatures,  and  back  pressures.  Proper  downstream  boundary  conditions  for  the 
subsonic  portion  of  the  outflow  are  shown  to  allow  variations  of  the  boundary  layer 
thickness  at  the  exit  plane  and  recirculating  separated  flows  for  sufficiently  high 
back  pressure.  Excellent  global  mass  conservations  are  demonstrated  by  using  the 
fully  conservative  form,  while  quasi-conservative  formulations  lead  to  unarrepi ably 
large  mass  conservation  errors. 

Along  with  the  investigations  of  Navier-Svckes  algorithms.'parabolized  Navier- 
Stokes(PNS)  procedures  are  also  studied.  The  PNS  al,:c  ihms  are  devised  from 
generalized  flux  split  TLNS  equations  which  include  both  the  traditional  pressure 
gradient  split  procedure  and  a  characteristics  split  system.  Comparisons  with  TLNS 


IV 


results  show  that  the  characteristica-baaed  PNS  system  gives  results  that  are  as 
accurate  as’  pressure-gradient-split  PNS  procedures.  The  use  of  a  safety  factor  in 
the  pressure  gradient  splitting  is  shown  to  cause  inaccuracies  and  should  be  avoided. 

The  global  pressure  iteration  for  the  PNS  algorithm  is  interpreted  as  an 
alternating-direction  procedure  for  the  TLNS  equations.  This  global  procedure 
is  shown  to  be  mathematically  well-posed  and  numerically  efficient. 

Swirling  viscous  flows  in  transonic  and  supersonic  propulsive  nozzles  have  been 
investigated  numerically,  The  central-difference  ADI  and  the  flux-vector  split  algo¬ 
rithms  are  utilized  to  solve  the  thin-layer  Navier-Stokes  equations  for  axisymmetric 
two-dimensional  flow  with  swirl.  The  effects  of  swirl  on  viscous  flow  are  studied  for 
nozzles  with  mild  to  high  expansion  ratios.  Both  flowfleld  detail  and  integral  nozzle 
performance  are  compared  to  previously  published  inviscid  calculations.  The  results 
show  that  the  presence  of  swirl  has  a  significant  effect  on  the  flowfleld  and  integral 
nozzle  performance,  especially  for  plug  nozzle  and  high  expansion  ratio  nozzles. 

Finally,  the  algorithms  developed  for  axisymmetric  two-dimensional  flows  are 
extended  to  solve  the  three-dimensional  TLNS  equations.  The  algorithms  are  based 
upon  DDADI  splitting  for  the  streamwise  flux  vector  and  additional  approximate 
factorization  of  the  operators  on  the  cross-stream  plane.  The  optimum  CFL  number 
reduces  to  the  order  of  10  and  it  gives  slower  convergence  as  compared  tr  corre¬ 
sponding  two-dimensional  algorithms  due  to  the  approximate  factorization  error. 
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CHAPTER  1 


INTRODUCTION 


Recent  interest  in  the  aerospace  plane  and  hypersonic  vehicles  has  revitalized 
research  on  high-speed  propulsion  systems.  In  the  design  of  a  propulsion  system, 
accurate  prediction  of  viscous  supersonic  flowfields  together  with  certain  physical 
parameters  such  as  thrust  play  a  critical  role.  Traditionally,  these  parameters  are 
obtained  from  wind  tunnel  tests  or  simplified  analytical  models.  The  analytical 
approach  is  only  valid  for  very  simple  geometries  and  flow  conditions  due  to  the  dif¬ 
ficulties  in  obtaining  exact  solutions  of  the  complicated  governing  equations.  Conse¬ 
quently,  successful  design  has  been  reliant  upon  expensive  wind  tunnel  experiments. 
With  the  advancement  in  computational  fluid  dynamics  (CFD)  and  computer  ar¬ 
chitectures,  numerical  computations  now  can  be  used  as  alternatives  of  experiments 
for  much  of  the  configuration  design  process.  Although  wind  tunnel  tests  continue 
to  be  important,  the  trend  is  clearly  toward  the  computational  approach  using  accu¬ 
rate  numerical  schemes  to  enhance  the  experimental  findings.  The  focus  of  present 
research  is  to  develop  accurate  numerical  algorithms  for  predicting  viscous  super¬ 
sonic  flowfields  that  occur  in  propulsion  systems.  In  particular,  the  predictions  of 
supersonic  flows  through  high  expansion  ratio  nozzles  will  be  emphasized. 

The  analysis  of  viscous  supersonic  flows  would  require  the  solution  of  the  com¬ 
pressible  Navier-Stokes  equations  with  proper  boundary  conditions.  It  is  well  known 
that  the  compressible  Navier-Stokes  equations  are  very  difficult  to  solve  because 
the  whole  equation  set  is  strongly  coupled  and  highly  non-linear.  To  avoid  directly 
solving  this  stiff  non-linear  system,  certain  degrees  of  approximations  have  to  be 
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made.  One  example  for  this  is  the  classical  Prandtl  boundary  layer  approach.  The 
boundary  Ihyer  assumptions  allow  inviscid  and  viscous  flows  to  be  computed  in¬ 
dependently.  For  supersonic  flows,  the  governing  equations  for  the  inviscid  region 
are  rendered  hyperbolic  by  neglecting  viscous  effects.  This  hyperbolic  equation  set 
can  be  efficiently  solved  by  the  method  of  characteristics  (MOC)  [lj.  For  the  vis¬ 
cous  region,  the  pressure  gradient  normal  to  the  wall  is  neglected  from  order  of 
magnitude  considerations;  thus,  the  Navier-Stokes  equations  reduce  to  boundary 
layer  equations.  Numerous  attempts  have  been  made  to  solve  the  boundary  layer 
equations  both  analytically  and  numerically.  The  analytical  technique  given  by  von 
Karman  and  Pohlhausen  (2]  requires  assumptions  of  the  velocity  profile  inside  the 
boundary  layer  and  is  only  valid  for  very  simple  problems.  The  numerical  solutions 
of  boundary  layer  equations,  which  can  handle  more  complex  problems,  have  been 
extensively  investigated  since  the  early  1970’s.  Some  representative  algorithms  are 
summarized  by  Anderson  et  at.  |3j. 

The  classical  boundary  layer  approach  assumes  the  interaction  between  the 
inviscid  region  and  the  viscous  region  is  small;  consequently,  either  region  can  be 
solved  independently.  To  take  into  account  this  interaction,  some  sort  of  inviscid- 
viscous  patching  procedure  has  to  be  employed.  One  typical  example  of  this  ap¬ 
proach  is  given  by  Ref.  (4).  The  patching  method  is  based  on  the  combination 
of  an  inviscid  MOC  procedure  and  a  boundary  layer  algorithm.  An  iterative  pro¬ 
cedure  between  inviscid  and  viscous  regions  is  accomplished  by  interchanging  the 
wall  pressure  from  the  MOC  procedure  and  the  displacement  thickness  from  the 
boundary  layer  procedure  until  convergence  is  achieved.  This  inviscid-viscous  inter¬ 
action  technique  does  pi  ’de  an  efficient  algorithm  to  calculate  viscous  supersonic 
flows.  However,  it  is  only  valid  for  weak-interaction  flows.  For  flows  with  strong 
interaction,  the  pressure  gradient  normal  to  the  wall  cannot  be  neglected.  Thus, 
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the  boundary  layer  algorithm  as  well  as  the  inviscid-viscous  interaction  procedure 
are  no  longer  applicable.  One  typical  example  of  these  strong  interaction  flows 
is  the  supersonic  flow  through  a  high-expansion  nozzle.  A  recent  work  done  by 
Kushida  (5j  indicates  that  the  boundary  layer  displacement  thickness  inside  the 
nozzle  can  be  as  large  as  42%  of  the  nozzle  radius  at  the  exit.  In  this  regime,  the 
inviscid-viscous  patching  procedure  fails  to  describe  the  pressure  variation  inside 
the  boundary  layer  and  the  realistic  mass  flow  rate,  thus  numerical  solutions  of  the 

\  v 

Navier-Stokes  equations  aje  required. 

For  typical  viscous  supersonic  flowfields,  the  governing  equation  set  is  hyper¬ 
bolic  in  the  supersonic  region  and  is  elliptic  in  the  subsonic  region  inside  the  bound¬ 
ary  layer.  This  mixed  hyperbolic/elliptic  character  makes  the  steady  Navier-Stokes 
equal  ions  cvtremcly  difF.CuIt  to  solve  because  a  different  numerical  algorithm  has  to 
be  trnpL/^d  hi  er.ch  region,  as  we  have  seen  in  the  classical  approach.  However,  if 
we  coi.s'der  the  unsteady  Navier-Stokes  equations,  the  equation  set  becomes  hyper¬ 
bolic  in  time  for  both  supersonic  and  subsonic  regions.  Therefore,  given  an  initial 
guess  of  the  ftowfield,  the  solutions  can  be  obtained  by  marching  in  time  until  the 
steady  state  is  reached.  This  procedure,  generally  referred  to  as  a  time-dependent 
or  time-iterative  scheme,  enables  one  numerical  algorithm  to  be  used  throughout 
the  flowfield.  The  time-dependent  concept  was  first  applied  to  inviscid  calcula¬ 
tions  for  flows  over  blunt  bodies  by  Mnretti  and  Abbett  |6|  in  19C6.  Since  then, 
time-dependent  solutions  have  become  an  important  segment  of  CFD.  The  first  ap¬ 
plication  in  compressible  viscous  flows  was  done  by  MacCormack  [7]  in  19C9.  In 
this  early  work,  an  explicit  predictor-corrector  scheme  was  proposed  to  solve  the 
Navier-Stokes  equations.  This  method  is  very  straightforward  to  program  but  it 
suffers  from  a  limitation  on  the  time  step  size  when  only  steady-state  solutions  are 
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Time-dependent  schemes  can  also  be  implemented  in  an  implicit  fashion.  The 
implicit  time-dependent  formulation  imposes  no  stability  limitation  on  the  size  of 
time  steps,  hence,  in  most  cases,  is  superior  to  the  explicit  scheme  if  only  steady- 
state  solutions  are  concerned.  One  important  application  of  implicit  time-dependent 
algorithms  to  compressible  Navier-Stokes  equations  is  the  alternating  direction  im¬ 
plicit  (ADI)  scheme  suggested  by  Beam  and  Warming  (8,9j,  which  is  also  noted 
as  the  linearized  block  implicit  (LBI)  scheme  by  Briley  and  McDonald  (10).  The 
ADI  scheme  has  gained  popularity  since  the  mid-70's  due  to  its  capability  to  solve 
multi-dimensional  inviscid  as  well  as  viscous  flows. 

W:th  the  progress  in  CFD  during  the  past  decade,  numerous  well-developed 
algorithms  are  now  available  for  compressible  Navier-Stokes  calculations.  These 
algorithms  an  in  general  be  divided  into  two  categories  according  to  the  the  type 
of  spatial  discretizations.  For  those  of  central-difference  type,  Steger  (11)  formulated 
the  ADI  scheme  in  the  general  coordinate  system,  Baldwin  and  Lomax  (12)  solved 
the  thin-layer  Navier-Stokes  (TLNS)  equations  with  an  algebraic  turbulence  model, 
and  Pulliam  (13)  applied  the  implicit  ADI  scheme  to  solve  flows  over  airfoils.  For 
those  of  upwind-difference  type,  Lombard  et  al.  ( 1 4 1  proposed  a  conservative  supra- 
characteristics  method  (CSCM)  based  on  non-conservative  flux-difference  splitting, 
and  MacCormack  developed  a  line  Gauss-Seidel  procedure  baaed  on  Steger  and 
Warming  ( 16)  flux-vector  splitting.  Similar  investigations  are  also  noted  by  other 
authors;  these  include  the  relaxation  scheme  by  Chakravarlhy  (17),  the  LU  scheme 
by  Yoon  and  Jameson  (18),  the  single  level  scheme  by  Lombard  et  al.  jl9|,  and 
the  diagonally  dominant  ADI  scheme  by  Chang  et  al.  (20|.  Thomas  and  Walters 
(2 1 1  used  a  similar  relaxation  procedure  to  solve  two-dimensional  viscous  supersonic 
flows  based  on  van  Leer’s  flux-vector  splitting  (22) .  More  recently,  this  work  has 
been  extended  to  three  dimensions  by  Newsome  et  al.  (23). 
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In  developing  an  efficient  numerical  algorithm  that  is  well  suited  to  the  vis¬ 
cous  supersonic  calculations  required  for  this  study,  four  aspects  of  solutions  of  the 
Navier-Stokes  equations  are  considered: 

1.  The  algorithm  should  be  abl  .o  take  into  account  the  predominantly  supersonic 
nature  of  the  flowfield,  and  consequently  give  rapid  convergence  in  the  high 
Reynolds  number,  unseparated  limit. 

2.  For  lower  Reynolds  number  flows,  proper  downstream  boundary  conditions 
have  to  be  implemented  on  the  subsonic  portion  of  the  exit  profile  so  that  the 
flow  will  respond  to  downstream  environmental  changes. 

3.  To  predict  thrust  with  accuracy,  global  mass  conservation  has  to  be  ensured. 
This  feature  is  accomplished  by  using  the  strong  conservative  form  of  the  gov¬ 
erning  equations. 

4.  The  algorithm  can  be  easily  simplified  to  a  certain  extent  such  that  a  pure 
space-marching  procedure  is  allowed  for  high  Reynolds  number,  "nseparated 
flows.  For  this  reason,  the  parabolized  Navier-Stokes  (PNS)  procedure  is  also 
considered  in  this  study. 

To  begin  with,  the  implicit  time-dependent  scheme  is  first  applied  to  solve 
the  quasi  one-dimensional  Euler  equations  for  spatial  discretizations  based  on  both 
central  differencing  and  upwind  differencing.  This  preliminary  work  allows  the  first 
assessment  of  algorithms  in  terms  of  computational  efficiency  and  accuracy.  Some 
details  of  the  algorithms  such  as  effects  of  approximate  Jacobians,  and  comparisons 
of  accuracy  between  first  order  and  second  order  upwind  schemes,  will  be  discussed. 

Two-dimensional  calculations  start  with  the  application  of  the  ADI  scheme 
to  the  axisymmetric  two-dimensional  TLNS  equations  in  order  to  justify  the  ap¬ 
propriateness  of  this  algorithm  for  viscous  supersonic  computations.  To  encounter 
predominantly  supersonic  flows,  a  hybrid  upwind/central  differencing  scheme  is  pro- 
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posed  along  with  its  Fourier  stability  analysis  (24).  Accordingly,  three  approximate 
factorization  algorithms  and  one  direct  method  are  formulated  for  the  solutions  of 
the  discretized  TLNS  equations  based  on  this  hybrid  differencing  scheme.  To  verify 
the  accuracy  of  the  proposed  hybrid  scheme,  the  results  computed  by  using  current 
algorithms  will  be  compared  to  those  by  the  MOC  procedure  given  in  Ref.  |-t]. 

For  better  understanding  of  the  effects  of  downstream  boundary  conditions 
on  the  flowfield,  supersonic  flows  through  a  conical  nozzle  and  a  high  area  ratio 
contoured  nozzle  are  computed  by  using  the  proposed  algorithms.  The  variation  of 
flow  character  is  obtained  by  varying  the  back  pressure  level.  In  particular,  back 
pressure  levels  that  are  sufficiently  high  to  produce  separation  i»  side  the  nozzle  are 
considered  in  order  to  simulate  the  classical  experimental  characteristics  that  are 
observed  when  altitude  nozzles  are  operated  on  sea-level  thrus'„  stands.  The  flowfield 
demonstrations  include  both  laminar  and  turbulent  calculations.  The  turbulent 
calculations  are  based  on  the  Baldwin  and  Lomax  model  1 12,25].  Comparisons 
of  global  mass  conservation  between  strong  conservative  and  weak  conservative 
formulations  are  made. 

Parallel  to  the  development  of  Navier-Stokes  algorithms,  the  applications  of  the 
time-dependent  scheme  on  PNS  procedures  are  also  studied.  Parabolized  Navier- 
Stokes  algorithms  |26-29]  have  proven  to  be  very  popular  because  of  their  accu¬ 
racy  and  efficiency.  For  many  flowfields,  they  give  results  that  are  almost  identical 
to  those  obtained  with  full  Navier-Stokes  equations,  although  the  CPU  lime  re¬ 
quired  is  much  less  than  that  needed  for  the  complete  equations.  The  basic  idea 
of  PNS  schemes  is  to  render  the  steady  state  Navier-Stokes  equations  parabolic  in 
the  streamwise  direction  by  proper  approximations.  This  parabolic  set  of  equations 
can  then  be  solved  by  a  space-marching  procedure  similar  to  the  MOC  procedure 
used  for  inviscid  supersonic  flows.  The  PNS  algorithms  differ  from  the  classical 
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boundary  layer  approach  in  that  the  normal  pressure  gradient  inside  the  boundary 
layer  is  retained  and  coupled  to  the  pressure  variation  of  the  inviscid  core  flow  in 
the  parabolized  equations.  Consequently,  PNS  schemes  can  handle  strong  inviscid- 
viscous  interaction  flows  without  losing  accuracy.  The  drawback  of  PNS  algorithms 
is  that  the  marching  procedure  fails  if  reverse  flow  is  present  in  the  flowfield  [3|. 

The  major  difference  between  PNS  procedures  and  Navier-Stokes  solvers  is 
that  PNS  schemes  are  normally  formulated  in  terms  of  the  steady  state  equations 
(see,  for  example  Ref.  (26))  while  Navier-Stokes  schemes  are  generally  formulated 
in  terms  of  the  time-dependent  equations.  Because  of  this,  it  is  difficult  to  extend  a 
PNS  algorithm  to  a  Navier-Stokes  algorithm.  In  the  present  study,  PNS  algorithms 
are  obtained  as  a  simplification  of  the  time-dependent  general  flux  split  Navier- 
Stokes  algorithms.  One  advantage  of  this  is  that  a  number  of  PNS  approximations 
can  be  defined  including  the  traditional  Vigneron  approach  (26)  and  a  new  approach 
based  upon  the  physical  characteristics  of  the  equations.  Furthermore,  the  resulting 
PNS  procedure  still  contains  the  temporal  derivative.  This  requires  the  solutions  to 
be  obtained  by  iterations  in  time  at  every  streamwise  station.  This  time-iterative 
PNS  procedure  makes  th«.  space-marching  problem  well-posed  and  consequently 
eliminates  the  necessity  of  a  safety  factor  that  occurs  in  the  *  aditional  approach. 

As  a  further  example  of  the  application  of  Navier-Stokes  solvers  mentioned 
above,  axisymmetric  swirling  nozzle  flows  are  studied.  Swirling  flows  ahead  of 
the  combustor  in  ramjet  applications  have  been  suggested  as  a  means  to  reduce  the 
reattachment  length  of  the  combustor  flowfield.  The  introduction  of  swirl  generated 
by  fixed  vanes  located  in  the  inlet  of  the  dump  combustor  can  greatly  increase  the 
efficiency  of  the  combustion  process  and  thus  reduce  the  length  of  the  combustor 
|30).  However,  the  residual  swirling  flow  in  the  combustor  will  enter  the  exhaust 
nozzle,  resulting  in  losses  in  thrust  and  reducing  the  mass  flow  rate.  Doth  of  these 
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decrease  the  nozzle  performance.  Therefore,  it  is  important  to  understand  to  what 
degree  swirling  affects  the  nozzle  flowfield  and,  subsequently,  the  overall  nozzle 
performance.  Several  previous  investigations  have  considered  the  effects  of  swirl, 
but  have  ignored  the  effects  of  viscosity.  In  this  study,  we  look  at  the  effects  of  swirl 
as  a  function  of  nozzle  Reynolds  numbers. 

Previous  investigations  of  swirling  nozzle  flow  include  both  quasi-one- 
dimensional  and  axisymmetpc  two-dimensional  analyses.  Carpenter  et.  al.  |3l) 
in  an  early  study  obtained  one-dimensional  results  by  neglecting  the  radial  veloc¬ 
ity  component.  Hoffman  and  co-workers  (32,33]  studied  swirling  flows  in  annular 
propulsive  nozzles  by  means  of  two-dimensional  inviscid  numerical  technioues.  To 
parameterize  their  studies,  they  used  four  different  inlet  swirl  profiles:  free  vortex, 
constant  angle,  forced  vortex,  and  Rankine  vortex.  Their  calculations  are  based 
upon  the  explicit  MacCormack  scheme  [7]  for  the  transonic  flowfield,  while  the 
method  of  characteristics  was  used  to  compute  the  supersonic  flowfield  after  the 
throat.  They  concluded  that  for  values  of  swirl  often  encountered  in  ramjet  and 
turbojet  applications,  the  effect  cf  swirl  on  the  nozzle  performance  is  small  and  ran 
probably  be  neglected. 

A  recent  work  by  Dutton  J3*i ]  indicates  that  significant  reductions  in  the  nozzle 
discharge  coefficient  and  the  vacuum  stream  thrust  efficiency  may  occur  for  hi^ii  val¬ 
ues  of  swirl  at  the  inlet  of  the  nozzle.  Again,  Dutton  usss  the  -  MacCormack 

scheme  to  analyze  three  different  nozzles,  including  a  convergent-divergent  (C-D) 
nozzle,  an  annular  nozzle,  and  a  converging  nozzle  Several  inlet  swirl  profiles  were 
enforced  as  inlet  boundary  conditions,  and  the  corresponding  effects  of  them  were 
identified.  He  also  verified  the  numerical  results  by  comparing  the  computed  wall 
static  pressure  with  experiments  for  a  C-D  nozzle  with  an  area  ratio  of  .25. 


9 


The  swirling  flow  investigations  mentioned  above  are  all  confined  to  inviscid 
calculations.  As  indicated  before,  the  boundary  layer  displacement  thickness  inside 
high  area-ratio  nozzles  can  be  very  large  at  the  exit.  In  this  regime,  the  inviscid 
assumption  is  inadequate.  The  present  study  proceeds  with  the  numerical  solutions 
of  viscous  swirling  nozzle  flows  by  using  implicit  time-dependent  schemes.  Viscous 
calculations  are  done  for  a  series  of  Reynolds  numbers  to  identify  the  effect  of  the 
boundary  layer  on  swirling  nozzle  flows.  To  place  these  viscous  results  in  perspective 
with  inviscid  calculations  appearing  in  the  literature,  the  results  in  the  inviscid 
limit  are  also  presented  along  with  those  of  tlv*  viscous  calculations.  Additional 
calculations  of  swirling  flows  in  high  expansion  nozzles  are  also  given.  Both  flowfield 
details  and  the  effect  of  swirl  on  the  integral  nozzle  performance  are  shown. 

Finally,  numerical  algorithms  developed  for  axisymmetric  two  dimensional 
flows  are  extended  for  three-dimensional  viscous  supersonic  calculations.  Both  PNS 
and  global  Navier-Stokes  procedures  are  demonstrated  by  flowfield  predictions  on 
a  three-dimensional  nozzle  with  a  rectangular  cross-section. 


■ 


CHAPTER  2 


THE  APPLICATION  OF  TIME-ITERATIVE  SCHEMES  TO 
THE  ONE-DIMENSIONAL  EULER  EQUATIONS 


This  research  starts  with  quasi  one-dimensional  calculations  of  compressible 
flows  for  two  reasons.  First,  the  analytical  solutions  of  these  flows  are  easily  ob¬ 
tained,  and  thus  provide  back-to-back  checks  of  the  accuracy  of  the  numerical  algo¬ 
rithms.  Second,  the  simplicity  in  the  formulation  of  the  equations  allows  a  series  of 
numerical  experiments  to  be  done  in  order  to  explore  the  potential  difficulties  asso¬ 
ciated  with  multi-dimensional  calculations.  The  central-differencing  as  well  as  the 
upwind-differenced  Euler  implicit  schemes  are  applied  to  the  calculation  of  quasi- 
one-dimensional  flows  through  a  convergent-divergent  nozzle.  Special  emphases  are 
placed  on  stability  analyses  of  the  numerical  algorithms  and  the  distinctive  charac¬ 
teristics  of  supersonic  flows. 

2.x  Gfly&rjiwg  gqygtippg 

The  unsteady  quasi  one-dimensiona!  Euler  equations  are  given  by 

|(,a)  +  £(,ua)=0 

± + a fH  =0  (2.1) 

d  d 

+  jj^\ (t  +  P)ual  =  0 

where,  standard  fluid  dynamic  notations  have  been  used.  These  include  the  density 
p ,  velocity  u,  pressure  p,  and  the  cross-sectional  area  a,  The  total  energy  e  per  unit 
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volume  is  defined  by 


1  i 

e  =  Pc  +  ~pu 


in  which,  c  is  the  internal  energy  per  unit  mass.  For  compressible  flow,  the  perfect 
gas  relation  is  used  to  close  the  problem. 

For  easier  implementation  of  numerical  procedures,  Eq.  (2.1)  is  expressed  in 
vector  form  as 

£  +  ^  =  *  <2'2> 
at  ax 

where  Q,  E ,  and  H  are  flow  variables,  flux  vector,  and  source  vector,  respectively. 
They  are  defined  by 

Q  =  \pa,pua,ta\ 

E  =  \pua,  ( pu 2  +  p)a,  (e  +  p)ua)T 

H  =  |0,,>j|,0| 

where  the  superscript  T  refers  to  the  transpose  of  the  vector.  Equation  (2.2)  is  writ¬ 
ten  in  strong  conservative  form  (35),  which  is  preferred  for  numerical  computations 
because  it  conserves  mass,  momentum,  and  energy  identically  in  the  discretized 
form.  For  flows  with  discontinuities,  this  conservative  formulation  allows  the  exis¬ 
tence  of  weak  solutions,  thus  allowing  shock-capturing. 

The  unsteady  Euler  equations  are  hyperbolic  in  time  and  can  be  converted  into 
uncoupled  characteristic  equations.  If  we  define  the  Jacobian  matrix  A  by 

dQ' 


and  use  the  chain  rule,  Eq.  (2.2)  becomes 

dQ  OQ 

+  A  ■—  =  H. 
dt  dx 

For  the  present  one-dimensional  case,  A  is  found  to  be 


(2.3) 


A  = 


0 

f-V 


+  i)« 


3  21  _  2 


1 

(3  -  -r)u 
-  i)«: 


o 

i  -  i 
*yu 
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The  Jacobian  matrix  A  can  be  transformed  to  a  diagonal  matrix  via  the  similarity 
transformation  defined  by 

A  =  MKM~ 1 .  (2.4) 

The  diagonal  matrix  A  takes  the  form 

A,  0  0  ‘ 

A  =  0  A3  0 

0  0  A3 

where  At,A3,and  A3  are  eigenvalues  of  the  matrix  A.  Matrices  M  and  M~l  are 
composed  of  the  left  and  right  eigenvectors  of  the  matrix  A,  respectively.  For  the 
matrix  A  given  above,  three  eigenvalues  are 

Aj  =  u 
Aj  =  u  +  e 
A3  =  v  -  c 

in  which,  c  is  the  speed  of  sound.  The  left  and  right  eigenmatrices  M  and  M" 1  are 
given  by 

’ 1  iu  iu 

M  =  U  ^j(c  +  !)  ^(c  “  !) 

.V  +  ^(37  +  -  U). 

and 

1  _  Lljlpl  L=p  ‘ 

+  o -nr)?)  ^  . 

l^<:Lr^+“>  ^;(-1  + (*-•>)?) 

Equation  (2.3)  now  becomes 


~  +  MAA/*'~ 
at  dx 


If  we  define  the  characteristic  variable  Q  by 


oq 

dt 
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and  multiply  Eq.  (2.2)  by  M  ,  we  have 


at  dx 


(2.S) 


where 


H  =  M~lH. 


Equation  (2.5)  is  equivalent  to  the  three  decoupled  characteristic  equations 

.  =  1,2,3  (2.6) 


dq,  .  dq,  * 


•  *  • 

in  which  and  h,  are  elements  of  Q  and  H,  respectively. 

The  procedure  above  demonstrates  that  the  one-dimensional  Euler  equations 
can  be  transformed  into  three  characteristic  equations  with  each  equation  governing 
one-dimensional  wave  propagation  with  a  specific  direction.  These  characteristic 
equations  can  be  obtained  by  multiplying  the  governing  equations  in  vector  form 
by  the  eigenmatrix  M~l.  For  subsonic  flows,  At  and  A2  are  positive  .while  As  is 
negative.  The  equation  set  possesses  both  right  and  left  running  characteristics. 
For  supersonic  flows,  all  three  eigenvalues  are  positive,  thus  the  waves  can  only 
travel  from  upstream  to  downstream.  As  will  be  discussed  later  in  this  chapter,  this 
allows  a  marching  procedure  to  be  used  for  supersonic  flow  calculations. 

2.2  The  Central-Differencing  Algorithm 

To  solve  Eq.  (2.2)  numerically,  the  central-differencing  Euler  implicit  scheme 
is  considered.  Symbolically,  the  Euler  implicit  scheme  can  be  expressed  as 


r»+  1 


Qn+l  _  q n  gg 

+  5=0 


(2.7| 


where  superscripts  n  +  1  imply  these  quantities  are  to  be  evaluated  at  the  new  time 
level.  If  we  define  AQ  =  Qn+1  -  Qn,  the  flux  vector  E  and  the  source  vector  II 
can  be  linearized  according  to  the  following  local  Taylor  series  expansions 

En  +  {  =  En  +  AAQ  (2.S) 
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//"+»  =  Hn  +  DAQ  (2.9) 

i 

in  which  D  it  the  Jacobian  matrix  defined  by  D  =  dH/dQ.  For  the  present  quasi 
one-dimensional  cate,  the  matrix  D  is  given  by 


D  = 


Upon  substitution  of  Eq.  (2.8)and  Eq.  (2.9)  into  Eq.  (2.7),  we  have 


(/  -  AtD  +  Al  —  A)A<?  =  -A tR 
ox 


(2.10) 


where  R  is  the  residual  vector  evaluated  at  time  level  n, 


(2.11) 


Note  that  all  the  derivatives  did i  in  Eq.  (2.10)  and  Eq.  (2.11)  imply  discretizations 
by  central  differencing.  The  left-hand  side  operator  of  Eq.  (2.10)  results  in  a  block 
tri-diagonal  matrix.  Each  block  is  a  3  x  3  matrix. 

2.2.1  Boundary  Conditions 

For  hyperbolic  equations,  the  boundary  conditions  can  be  easily  enforced  by 
using  the  MOC  boundary  procedure  suggested  by  Rai  and  Chaussee  130  and 
Chakravarthy  (37) .  As  indicated  earlier,  the  governing  equations  imply  three  waves 
travel  with  specific  directions.  Boundary  conditions  are  imposed  for  those  waves 
running  into  the  computational  domain,  while  for  waves  moving  from  inside  (lie 
domain  toward  the  boundary,  the  decoupled  characteristic  equations  snrli  a*  those 
given  in  Eq.  (2.6)  are  used  to  allow  the  information  to  propagate  from  inside  the 
domain. 

For  subsonic  inflows,  A !  and  A  j  are  positive,  which  implies  two  conditions  must 
be  specified  at  the  upstream  end.  A  reasonable  choice  is  to  specify  the  stagnation 


pressure  P°  and  the  stagnation  temperature  T°.  Let  these  specified  values  of  P° 
and  T°  be  given  as 


P°  =  K\ 
T°  =  K7. 


If  we  define  a  vector  fi  by 

n  =  (p°,r°,o)T, 


then,  from  Taylor  series  expansion  of  fi,  we  ha-'e 

Afi  =  nn+1  -  nn  =  |^aq 

where  dU/dQ  is  the  Jacobian  matrix  of  fi.  To  force  fi'*+1  to  be  fixed  at  the  value 
of  fi,,  where  fi,  =  (ft|,Kj,0)T,  the  following  equation  can  be  employed 

^AQ  =  n,-nn,  (2.12) 


Since  the  third  eigenvalue  A3  is  negative,  we  must  select  the  decoupled  char¬ 
acteristic  equation  corresponding  to  A3  from  Eq.  (2.6)  to  complete  the  upstream 
boundary  conditions.  If  we  define  the  selection  matrix  L~  by 


L~  = 


0  0 

0  0 

0  0 

m 


o' 

0 

1 


and  multiply  L~M~l  on  both  sides  of  Eq.  (2.10),  we  obtain  the  characteristic 
equation  corresponding  to  A3  as 


M  '(/  -  A  ID  +  At—  /t)AQ  =  -btL'Xr'R. 
ox 


(2.13) 


Combining  Eq.  (212)  and  Eq.  (2.13),  the  discretized  equation  at  the  inlet  boundary 
can  be  written  as 


,<?fi  .. 

0Q  + 


A/“ '(/  -  AtD  +  At^/1))AQ  =  fi,  -  fi"  -  AtL~Ar '/?. 
ax 


(2.1-1) 
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For  supersonic  flows  at  the  inlet,  all  three  characteristics  come  from  outside 
the  domain;  therefore  all  entries  of  the  dependent  variable  Q  have  to  be  specified. 
For  supersonic  flows  at  the  exit,  all  three  characteristics  are  outgoing;  therefore,  the 
discretized  equation  itself,  Eq.  (2.10),  can  be  applied  directly  without  any  special 
treatment.  For  subsonic  flows  at  the  exit,  A3  is  negative,  hence  one  boundary 
condition  has  to  be  specified.  Let  the  specified  quantity  be  the  back  pressure  = 
K3.  The  vector  fi  now  takes  the  form 

n  =  (o,o,p6)r. 

To  select  the  characteristic  equations  corresponding  to  Aj  and  Aj,  one  can 
choose  the  selection  matrix  L+  as 


L+  = 


1  0 

0  1 

0  0 

» 


Similarly,  the  discretized  equation  for  this  case  is 

l|^  +  L+ATl(I  -  UD  +  At|^))AQ  =  n*  -  nn  -  AtL+Al'lR  (2.15) 

where  the  constant  vector  fle  =  (0,0,  K3)  . 

In  the  discretized  equations  at  the  boundaries,  Eq.  (2.14)  and  Eq.  (2.15),  the 
centrally  differenced  spatial  derivatives  d/dx  are  not  applicable.  To  remedy  this, 
we  use  two-point  one-sided  differences  instead  of  central-differences.  This  approach 
retains  the  block  tri-diagonal  structure  of  the  left  hand  side  matrix  but  is  only  first 
order  accurate. 

In  order  to  have  better  solution  accuracy,  three-point  one-sided  differences  can 
be  used.  This  results  in  extra  elements  at  the  first  and  the  last  row  of  the  left 
hand  side  matrix,  which  can  be  eliminated  easily  by  elementary  matrix  operations 
(38) .  This  approach  retains  second  order  spatial  accuracy  throughout  the  whole 
computational  domain,  and  will  be  generally  used  for  the  discretized  equations  at 
the  boundaries. 


17 


2.2.2  Stability  Analysis 

The  application  of  Fourier  or  von  Neumann  stability  analysis  (24)  has  become  a 
powerful  tool  for  today’s  CFD.  In  developing  a  new  numerical  algorithm,  the  stabil¬ 
ity  analysis  provides  abundant  information  about  the  convergence  requirements  of 
various  parameters  involved  in  the  algorithm.  Before  attempting  to  solve  Eq.  (2.2) 
numerically,  we  consider  the  Fourier  analysis  of  its  discretized  form,  Eq.  (2.10). 

For  any  given  function  /(x,t),  the  Fourier  transform  is  defined  by 


/(uM)  =  f{x1t)e~'u,tdx 

J  —  OO 


where  i  is  the  square  root  of  --1.  This  transform  exists  only  if  f(x,t )  is  square 


summable,  that  is, 


/OO 

-OO 


X  <  OO. 


The  inverse  transformation  which  transforms  /  from  the  frequency  domain  to  the 


spatial  main  is  defined  by 


/(*,<)  =  I' /(w.iK-jj. 


The  analogous  transform  for  a  function  q{x,t)  defined  only  at  discretized  points  can 


be  written  as 


r»  =  AX  £ 

»=  — OO 

where  the  superscripts  n  denote  the  time  step  ( t  =  nAt)  and  the  subscripts  i 
represent  the  spatial  step  (r,  =  i  Ax).  The  inverse  transformation  for  the  discretized 


function  q  is  given  by 


r/Ai  du 

1?  =  / 

J -n/dz 


(2.10) 
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Upon  substitution  of  a  specific  Fourier  mode  with  frequency  w  into  the  dis¬ 
cretized  equation,  we  obtain  the  functional  relationship  of  the  amplification  factor 
defined  by 


iM  1 


M 


9  -  W '  <217> 

The  stability  criteria  for  any  specific  algorithm  are  then  determined  by  the  magni¬ 
tude  of  g.  If  |p|  is  greater  than  unity,  the  amplitude  corresponding  to  the  wave  mode 
u  is  growing,  and  hence  is  unstable.  If  \g\  is  less  than  unity  for  all  wave  modes,  the 
algorithm  is  stable. 

In  solving  the  central-difference  discretized  equation  numerically,  second  order 
as  well  as  fourth  order  artificial  dissipation  terms  are  added  to  Eq.  (2.2)  to  avoid 
odd-even  decoupling  and  to  damp  out  high  frequency  oscillations.  This  results  in 


dQ  +  dJL  _  liA  2  d*Q  =  u  _  <«Ar4  d4Q 

dt  dx  4  d2xdt  8  At  dx 4 


where  c,  and  <*  are  positive  constants.  The  discretized  equation,  Eq.  (2.10),  now 
becomes 


(/  -  AtD  +  At  9  A  -  liAx3~-/)A<?  =  -AtR  -  (2.18) 

ox  4  dx2  8  dx* 

For  linear  stability  analysis  A  can  be  treated  as  a  constant  matrix.  In  the 
frequency  domain  Eq.  (2.18)  becomes 

LiQn+l  =L2Qn  (2.10) 


where  L\  and  L2  are  given  by 


Lx 


-  +  ~(l  ~  cosu»x)|  +  t— - —  ,4sinu/.T 

**  L\  X 


-  AtD 


L i  -  / 1 1  f  --(1  -cosu^)  -  -‘(cosu'j. 


and 
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with  u >t  representing  the  wave  number  defined  by  ujs  =  wAx. 

Analogous  to  the  definition  of  the  amplification  factor  for  the  scalar  system, 
we  can  define  the  amplification  matrix  G  by 

Qn+l  =  GQn. 


The  convergence  criteria  are  then  determined  by  the  eigenvalues  of  the  matrix  G. 
A  stable  algorithm  is  ensured  when  the  magnitudes  of  eigenvalues  of  G  are  all  less 
than  unity.  From  Eq.  (2.19),  G  can  be  easily  evaluated  by  G  =  Ll[l L^.  At  the  high 
wave  number  limit  (w£  -  .. the  eigenvalues  of  G  are  found  to  be 


<7i 

02 

03 


1  4- 1,  -  2ct 
1  +  <« 


1  +  t,  -  2 £e 

1  +  <«  - 
l  +  «, 


where  a  is  defined  by 


a  = 


uAt 

~Kz 


which  is  referred  to  as  the  Courant-Friedriches-Lewy  (CFL)  number.  According  to 
the  absolute  values  of  g i  and  <73,  it  is  required  that 


0  <  l'  <  1  +  £, 


(2.20) 


to  maintain  numerical  stability.  The  value  of  gi  depends  on  the  CFL  number,  < , , 
tt,  and  the  geometry.  The  stability  criteria  associated  with  it  are  rather  involved. 
However,  several  conclusions  still  can  be  drawn.  First,  in  a  divergent  portion  of  the 
geometry  (jjf  >0),  is  always  less  than  unity,  thus  the  Euler  implicit  scheme 
is  stable.  Second,  in  a  convergent  section  <  0),  there  exists  a  certain  range 
of  CFL  such  that  |<7j|  >  1  for  fixed  t,  and  te.  Third,  if  no  fourth  order  artificial 
dissipation  is  included  (ce  =  0),  (92 1  is  always  greater  than  unity  in  the  divergent 
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section,  hence  the  algorithm  is  unstable.  In  practical  situations,  the  geometry 
contains  both  convergent  and  divergent  sections;  therefore,  the  central-difTerence 
Euler  implicit  scheme  for  quasi  one-dimensional  flows  is  only  conditionally  stable. 

The  eigenvalues  of  G  can  be  calculated  numerically  for  various  wave  numbers. 
Figure  1  shows  the  plot  of  the  magnitude  of  the  maximum  eigenvalue  versus  ut 
for  CFL  =  1, 10, 100  at  a  flow  Mach  number  of  0.5.  This  figure  clearly  illustrates 
that  increasing  the  CFL  number  tends  to  decrease  the  magnitude  of  the  maximum 
eigenvalue,  which  is  beneficial  for  the  speed  of  convergence.  Effects  of  the  artificial 
dissipation  are  demonstrated  in  Fig.  2,  where  the  maximum  eigenvalue  of  G  is 
plotted  for  <f  =  0,0.25,0.5, 1.0  for  a  fixed  CFL  of  10.  It  shows  that  the  addition  of 
fourth  order  dissipation  damps  out  high  frequency  components  of  the  wave.  It  is 
also  observed  that  fe  =  0.5  is  optimal  as  far  as  convergence  is  concerned. 

The  Fourier  analysis  discussed  above  is  based  on  two  major  assumptions.  First, 
the  analysis  is  only  valid  for  linear  cases,  in  other  words,  the  nonlinear  effect  of  the 
Jacobian  matrix  A  has  been  neglected.  Second,  the  analysis  assumes  an  infinite 
domain  and  excludes  the  effect  of  boundary  conditions.  Therefore,  the  results  are 
qualitatively  rather  than  quantitatively  accurate. 

From  the  results  of  Fourier  stability  analysis,  it  is  apparent  that  the  CFL 
number  plays  an  important  role  on  the  speed  of  convergence.  To  obtain  optimum 
convergence,  the  CFL  number  should  be  as  large  as  possible  provided  that  numerical 
stability  is  retained.  Since  the  CFL  number  is  directly  related  to  the  time  step  size, 
we  can  calculate  At  according  to  the  desired  CFL  number.  If  one  is  interested  in 
accurate  solutions  during  a  transient,  At  must  be  uniform  throughout  the  flowfield. 
In  this  case,  At  is  better  determined  according  to  the  maximum  value  of  u,  that  is. 

crAz 

At  =  - . 

Umar 

In  general,  the  Velocity  u  varies  from  point  to  point,  the  using  of  a  uniform  time  step 
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will  result  in  non-uniform  CFL  numbers  throughout  the  flowfield.  Consequently, 
the  overall  convergence  is  deteriorated,  especially  when  only  steady  state  solutions 
are  of  interest.  In  order  to  have  optimum  convergence,  At  can  be  locally  determined 
by  the  given  CFL  number,  in  other  words,  the  time  step  size  at  each  grid  point  is 
calculated  according  to 

,  a  Az 

At  =  - . 

u 

This  implies  a  constant  CFL  number  has  been  enforced  over  the  whole  flowfield.  The 
introduction  of  this  spatially  varying  time  step  (or  so-called  constant  CFL^  greatly 
enhances  the  speed  of  convergence  [39].  in  this  study,  the  constant  CFL  approach 
will  in  general  be  used  for  all  calculations  since  only  steady  state  solutions  are 
concerned. 

2.2.3  Computational  Results 

The  one-dimensionai  flow  through  a  convergent-divergent  nozzle  with  the  area 
variation  given  by 

.  2ttz  4-  '“it/i 

2  XL  2 

is  chosen  as  a  test  problem  The  geometry  associated  with  definitions  of  x,  Axrx, 
Ath ,  and  L  are  shown  in  Fig.  3.  A  uniform  grid  with  total  of  40  points  and  an  area 
ratio  ( Ath/A,n )  of  0.8  are  used  for  all  calculations  that  follow. 

Three  typical  cases  are  investigated,  including  pure  subsonic,  transonic,  and 
pure  supersonic  flows.  Figure  4  shows  the  convergence  rates  of  these  cases  by 
plotting  the  L-2  norm  of  the  non-dimensional  change  >n  Q  (&Q/Q)  against  the 
number  of  iterations.  Each  curve  is  obtained  by  using  the  optimum  CTL  number 
(the  CFL  number  that  gives  the  fastest  convergence).  The  supersonic  case  converges 
the  fastest  among  the  three  cases  because  of  the  predominantly  hyperbolic  nature 
of  the  flowfield.  For  pure  subsonic  flows,  the  system  of  equations  is  elliptic  in 


a(x)  = 


Subsonic 
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the  spatial  direction  since  both  left-running  and  right-running  characteristics  are 
present  in  the  flowfield.  The  subsonic  case  shows  the  slowest  convergence  due  to  this 
elliptic  behavior.  The  accuracy  of  the  central-difference  formulation  is  demonstrated 
in  Fig.  5  and  Fig.  6,  where  the  computed  Mach  number  and  pressure  distributions 
along  the  streamwise  direction  for  the  transonic  calculation  are  compared  to  those 
from  exact  solutions  for  the  same  nozzle.  The  comparison  shows  that  computational 
results  agree  very  well  with  exact  solutions. 

2.3  The  Upwlnd«Diflerenclng  Algorithm 

The  main  purpose  of  this  study  is  to  develop  efficient  numerical  algorithms  for 
supersonic  calculations.  As  indicated  earlier,  the  predominantly  hyperbolic  nature 
of  supersonic  flows  distinguishes  themselves  from  transonic  and  subsonic  flows.  To 
take  advantage  of  this  character,  upwind  schemes  appear  to  be  attractive.  As  we 
have  seen  in  Section  2.1,  the  Jacobians  of  the  governing  equations  generally  contain 
both  positive  and  negative  eigenvalues.  The  sign  of  each  eigenvalue  implies  the 
direction  of  wave  propagation  on  the  x  -  t  plane.  The  crux  of  flux-vector  splitting 
upwind  algorithms  ( 16,22)  is  to  separate  the  flux  vector  E  into  parts  with  definite 
(positive  and  negative)  eigenvalues.  The  splitting  can  be  formally  indicated  as 

E  =  E+ +  E~  (2.21) 

where  the  eigenvalues  of  the  Jacobian  of  E+  are  positive  and  those  of  E~  are 
negative.  There  are  an  infinite  number  of  ways  to  accomplish  this  splitting.  As 
an  example,  we  have  considered  the  Steger  and  Warming  (16)  splitting.  From  the 
similarity  transformation  of  A  defined  by  Eq.  (2.4),  we  readily  have 

A  = 

The  Steger  and  Warming  splitting  takes  the  form 

A+  =  (A  +  |  A  |) /2 
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Figure  5,  Mach  number  distribution  of  transonic  1-D  calculation, 
central-difference  solutions 
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Pressure  distribution  of  transonic  1-D  calculation  , 
central-difference  solutions 


A”  =  (A  -  |A|)/2 


where  |A|  refer#  to  the  matrix  composed  of  the  absolute  values  of  the  elements  of 
A.  The  matrix  A  can  then  be  decomposed  into 

A+  =  A/A+A/-‘ 

A~  =  A/  A  “  M  ~ 1 

with  A  A+  +  A~ .  By  using  the  homogeneous  property  of  the  matrix  A ,  we  have 

E  =  (A+  +  A~)Q  =  E+  +  E~ 


in  which,  E+  —  A+Q  and  E  —  A~  Q. 

Again,  using  Euler  implicit  differencing  in  time,  the  flux  split  system  can  be 


described  by 


Qn+'-Qn  ,dE+  dE~  .„n+1  .. 

+  (  —  +  —  _  H)  =0. 


(2.22) 


A  linearization  similar  to  Eq.  (2.8)  can  be  applied  to  E +  and  E~ .  The  result  ing 


delta  form  is 


|7  -  A.'D  +  +  ^-A;)\£Q  =  -At/?' 

Ox  Ox 


(2.23) 


where  the  residual  vector  R'  is 


,  dE+  dE~  " 


(2.24) 


The  Jacobian  matrices  At+  and  AJ  are  defined  by  A?  =  dE+  jdQ  and  /l'  - 
dE~  /dQ.  Note  that  A,  A  +  and  Aj"  ^  A~  when  the  eigenvalues  of  A  are  of 
mixed  sign.  For  the  one-dimensional  case,  if  the  flow  is  supersonic  (u  >  r).  ,1/  is 
exactly  the  same  as  A  given  in  Section  2.1  and  Af  is  zero.  As  an  approximation, 
A  +  and  A~  can  be  used  instead  of  At+  and  A,  in  Eq.  (2.23).  The  effect  of  using 
true  oi  approximate  Jacobian  matrices  is  detrimental  as  will  be  discussed  in  the 
next  section. 
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All  spatial  derivatives  in  Eq.  (2.23)  and  Eq.  (2.24)  imply  that  they  will  be 
differentiated  according  to  the  signs  of  their  eigenvalues.  For  example,  dE+  /dx 
and  dE~ /dx  are  differentiated  according  to 

dE+  E+-S+,  ,  E*  -  2Etx  +  S+ , 

dx  Ax  K  2A* 

and 

dE~  _  Ej+\  ~  Ej  _  Ei  —2 E^¥lEtA.7 
dx  Ax  2  Ax 

where  k  =  0  for  the  first-order  scheme  and  k  =  1  for  the  second-order  scheme. 

The  left-hand  side  matrix  in  Eq.  (2.23)  is  block  tri-diagonal  if  all  spatial  deriva¬ 
tives  are  discretized  by  first-order  upwind  differences.  It  is  well  known  that  first- 
order  upwind  differencing  adds  a  large  amount  of  artificial  dissipation  to  the  nu¬ 
merical  algorithm,  and  is  highly  inaccurate.  If  second-order  accurate  differencing  is 
employed,  the  left-hand  side  matrix  becomes  block  penta-diagonal,  which  is  more 
time-consuming  to  solve  than  the  block  tridiagonal  matrix;  however,  if  only  steady- 
state  solutions  are  of  interest,  one  can  use  first  order  differencing  on  the  left-hand 
side  and  second-order  differencing  on  the  right-hand  side.  This  will  retain  the  block 
tri-diagonal  structure  of  the  left-hand  side  matrix  while  maintaining  the  second- 
order  accuracy  of  the  steady  state  solution.  Jespersen  and  Pulliam  |40|  have  shown 
that  this  non-consistent  first  and  second  order  differencing  will  reduce  the  stability 
bound  of  the  CFL  number  and  slow  down  the  convergence.  To  make  consistent 
second-order  differencing  possible,  one  can  approximately  factorize  the  left-hand 
side  of  Eq.  (2.23)  as, 

a  a 

(/  -  AtD  +  Al—  Af)[I  -  AlD)'l{I  -  AtD  +  At—-A7)AQ  =  -AtR'  (2.25, 
dx  dx 

Equation  (2.25)  is  equivalent  to 

a 

(/  -  AtD+  At  —  At)AQ’  ^  -AtR' 
dx 


(2.2G) 
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and 

(/-  MD  +  Af^-Af)A<?  =  (/ -  A<Z?)A<?‘.  (2.27) 

ox 

Equation  (2.26)  can  now  be  solved  by  space-marching  from  upstream  to  downstream 
since  the  left  hand  side  matrix  is  lower  bi-diagonal  for  first-order  spatial  differenc¬ 
ing  and  is  lower  tri-diagonal  for  second-order  spatial  differencing  .  This  forward 
marching  together  with  the  backward  marching  given  in  Eq.  (2.27)  will  complete 
one  iteration. 

Two  comments  can  be  made  at  this  point.  First,  as  indicated  by  Sieger  and 
Warming  {16],  small  oscillations  occur  by  using  this  splitting  when  a  sonic  line  is 
crossed  because  of  the  discontinuity  in  the  first  derivative  of  the  split  flux  when  the 
eigenvalues  change  sign.  This  oscillation  can  be  partly  removed  by  the  introduction 
of  a  blending  coefficient  in  the  calculation  of  eigenvalues  [41].  Second,  if  the  flow 
is  supersonic,  E~  and  A,  are  identically  tero.  Equation  (2.23)  reduces  to  a  form 
analogous  to  the  forward-sweep  step  given  in  Eq.  (2.26),  namely 

a  a  r+  n 

(1  -  Aio  +  bt—At)AQ  =  -AH-JJ-  -  H)  .  (2.28) 

As  in  Eq.  (2.26),  the  left  hand  side  matrix  of  Eq.  (2.28)  is  lower  triangular  and  only 
the  forward  marching  step  of  Eq.  (2.26)  is  necessary  to  complete  one  iteration.  This 
allows  a  pure  marching  solution  to  be  obtained.  In  fact,  we  can  rearrange  Eq.  (2.28) 
to 


\l-AtD  +  i.l  +  l)~A;\£>Q,  =  -A<(~ 


+  K 


K_  -  2 c,  4-  e;_2 

2Ax 


-//) 

(2.29) 


where  k  =  0  for  Prst-order  scheme  and  k  =  1  for  second-order  scheme.  This 
equation  is  a  marching  equation,  which  can  be  solved  by  iterating  in  time  at  each 
grid  point  before  advancing  to  the  next  streamwise  location.  Thus,  the  vector  &Q, 
on  the  left  hand  side  of  Eq.  (2.29)  can  be  driven  to  the  desired  tolerance  by  time¬ 
marching  at  the  i-th  grid  point  before  the  procedure  marches  to  the  i+l-th  point. 
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This  marching  procedure  will  be  extended  for  two-dimensional  calculations  in  the 
following  ch'apters  and  will  not  be  discussed  in  detail  here. 

2.3.x  Boundary  gandlilgai 

The  boundary  procedures  for  upwind  schemes  are  simitar  to  those  for  central- 
difference  schemes.  At  boundaries,  the  characteristics  coming  from  outside  the 
domain  are  not  defined  and  are  replaced  by  specified  boundary  conditions.  By 
neglecting  these  incoming  characteristics,  the  discretized  equations  reduce  to 

(/  -  AtD  +  At  ~~A~)AQ  =  -  At(^ - H)  (2.30) 

OX  ox 

at  the  upstream  inlet,  and 

a  a  n 

(/  -  A tD  +  At— A*)AQ  =  -At(  —  -  H)  (2.31) 

at  the  downstream  end.  In  these  two  equations,  A+  and  A~  have  been  chosen 
instead  of  At+  and  Af  to  enable  the  application  of  the  MOC  procedure.  By  us¬ 
ing  the  identity  matrix  /  =  MM~X  and  the  definition  of  A  +  ,  Eq.  (2.30)  can  be 
approximated  by 

<>  n-  W 

A/(7  -  AtD  +  At  — A")Af"1  A(?  =  -  At(-^ - H)  . 

OX  ox 

Multiplying  both  sides  by  M_l,  we  have 

(/-  AtD  +  At^-A“)A<?  =  - AtA/"*(~—  -  H) 

OX  OX 

which  is  equivalent  to  the  previously  defined  decoupled  characteristic  equation. 
Again,  multiplying  this  characteristic  equation  by  a  selection  matrix  L+  and  com¬ 
bining  it  with  the  specified  boundary  conditions  discussed  in  Section  2.2.1  gives 

+  L+A/-‘(/  -  AtD  +  At  -~A~)\AQ  =  fi,  -  H"  -  AtL+ M~l  -  H) 


(2.32) 
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where  the  vectors  fl  and  fli  follow  the  same  definitions  given  in  Section  2.2.1.  This 
boundary  procedure  is  similar  to  Eq.  (2.15)  except  different  decoupled  discretized 
equations  have  been  used  in  Eq.  (2.32).  The  discussion  above  illustrates  that  the 
boundary  procedures  formulated  for  the  central-differencing  algorithm  given  in  Sec¬ 
tion  2.2.1  are  also  applicable  to  the  upwind-differencing  scheme. 

2.3.2  st3bilUy  A.p.ah:gis 


The  stability  analysis  is  now  studied  for  the  upwind  algorithm.  For  demonstra¬ 


tion,  only  first-order  differencing  on  both  sides  of  Eq.  (2.23)  is  considered.  Following 


the  definitions  given  in  Section  2.2.2,  the  amplification  matrix  G  can  be  expressed 


by  G  =  L\x  Li  with 


At  ,  At 

Li  =  /  -  A tD  +  — —  (1  -  cos  wx  4-  iai.iut)A7  4-  —  (cosw*  +  i sinu>x  -  l)Af 
Ax  Ax 


Li  -  I 


where  true  Jacobian  matrices  A*  and  /if  have  been  used.  Alternatively,  if  the  ap¬ 


proximate  Jacobian  matrices  A+  and  A~  are  used  on  the  left  hand  side  of  Eq.  (2.23), 


L|  and  Li  become 


At 

L\  —  I  -  AtD  4-  — —  (1  -  cosu/j  +  isinwx)/l  + 
Ax 

At 

+  - — (cosw*  4-  isinwx  -  1)A- 
Ax 

Li  =  I  4-  ~(1  -  cos  ux  +  isinwx)(/4+  -  /if) 
At 

4-  —  (cosu/x  4-  «sinu/x  -  1)(A~  -  /if). 


As  noted  earlier,  A?  and  Af  differ  from  A+  and  A~  when  the  flow  is  subsonic. 
Typical  stability  results  for  the  approximate  Jacobian  case  are  shown  in  Fig.  7  for 
a  flow  Mach  number  of  0.5.  An  explicit-like  CFL  restriction  ( o  <  1)  is  observed. 
On  the  other  hand,  if  true  Jacobians  are  used  on  the  left  hand  side  of  Eq.  (2.23), 
the  upwind  algorithm  is  unconditionally  stable,  as  shown  in  Fig.  8.  For  supersonic 
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flows,  A  f  =  A  and  Aj"  is  identically  zero.  Corresponding  stability  results  are  shown 
in  Fig.  9  for  a  flow  Mach  number  of  2.0.  The  upwind  scheme  is  again  stable  for 
all  CfL  numbers.  As  can  be  seen  in  Fig.  8  and  Fig.  9,  for  a  given  CFL  number, 
the  max.mum  eigenvalue  of  G  reaches  its  minimum  value  at  the  high  wave  number 
limit  (u>x  —  it).  This  characteristic  implies  that  the  upwind  scheme  is  naturally 
dissipative  and  no  artificial  viscosity  is  necessary  to  maintain  numerical  stability. 

2.3.3  Computational  Results 

The  same  test  problem  given  in  Section  2.2.3  is  calculated  by  using  Eq.  (2.23). 
Again,  three  typical  cases  are  studied,  they  are  pure  subsonic,  transonic,  and  pu»*e 
supersonic  flows.  Figure  10  compares  convergence  rates  obtained  by  using  first- 
order  differencing  on  both  sides  of  Eq.  (2.23)  for  all  cases.  It  shows  that  the  upwind 
algorithm  is  equally  efficient  as  the  central-difference  algorithm  (compare  Fig.  4), 
except  for  the  transonic  case,  for  which  the  discontinuity  across  the  sonic  point 

i 

substantially  slows  down  the  conveigence.  The  very  slow  convergence  of  the  sub¬ 
sonic  case  based  on  the  approximate  Jacobians  where  the  optimum  CFL  number  is 
found  to  be  0.9  is  also  shown  on  Fig.  10.  Tn  fact,  with  the  use  of  the  approximate 
Jacobian,  the  computer  code  diverges  for  both  transonic  and  subsonic  calculations 
•f  o  >  1,  thus  confirming  the  stability  predictions  given  in  the  last  section. 

Computational  results  lor  the  supe/sonic  calculation  are  compared  to  exact 
solutions  in  Fig.  11  and  Fig.  12  for  both  first  order  and  second  order  accurate  com¬ 
putations.  The  second-order  scheme  gives  solutions  that  are  much  more  accurate 
than  the  first-order  scheme  doe-.  Therefore,  second-order  differencing  should  always 
be  used  to  ensure  accurate  solutions. 


Subsonic  approximate  Jacobian 


Number  of  iterations 


Figure  10.  Convergence  for  1-D  implicit  first-order  upwind  scheme 
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Figure  11.  Mach  number  distribution  of  supersonic  1-D  calculation, 
upwind  solutions 


CHAPTER  3 


THE  APPLICATION  OF  TIME-ITERATIVE  SCHEMES  TO 
THE  TWO-DIMENSIONAL  NAVIER-STOKES  EQUATIONS 


Based  on  the  information  gained  from  the  application  of  time-marching  schemes 
to  the  one-dimensional  Euler  equations,  this  chapter  proceeds  with  numerical  so¬ 
lutions  of  the  axisymmetric  two-dimensional  Navier-Stokes  (N-S)  equations.  The 
conventional  implicit  ADI  procedure  is  first  applied  to  transonic  and  supersonic 
viscous  calculations.  The  appropriateness  of  this  procedure  when  applied  to  pre¬ 
dominantly  supersonic  flows  wilt  be  identified.  According  to  the  physics  of  viscous 
supersonic  flows,  a  discretized  scheme  using  upwind  flux-vector  splitting  in  the 
streamwise  direction  and  central-differencing  in  the  cross-stream  direction  together 
with  solution  procedures  are  proposed.  The  Fourier  stability  analysis  will  be  used 
to  analyze  the  stability  criteria  of  this  new  discretized  scheme.  Of  the  solution 
procedures,  approximate  algorithms  as  well  as  a  direct  solver  are  considered.  These 
procedures  will  be  used  to  calculate  viscous  supersonic  flows  through  nozzles.  In 
particular,  attention  will  be  paid  to  proper  downstream  boundary  conditions  for 
the  subsonic  portion  of  the  outflow  and  global  mass  conservation. 

3.1  Governing  Equations 

For  practical  applications,  the  two-dimensional  Navier-Stokes  equations  are 
formulated  in  a  cylindrical  coordinate  system.  The  equations  for  planar  two- 
dimensional  flows  can  be  easily  obtained  by  simplifying  the  cylindrical  version. 
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Let  z  and  y  denote  the  axial  and  radial  coo  '  The  unsteady  axisymmel 


ric  Navier-Stokes  equations  for  laminar  <or  are  given  as 


d  ,  ^  0  ,  v  8  X  A 

gjM  +  ^  -  0 


f  (pop)  +  £(/»«»)  +  ^l(P"2  +  P)vl  =  P  +  •  V  -  +  §j)»l 


\ 

£'i‘v)  +  jjl(*  +  p)»»l +  5l(«  +  p)»»l  -  Tz (l"u,Jll  -  5 V ■  V)  +  «»(|j 


(3.1) 


du.  dT .  .  d  .  dv  du.  dv  2_  7.  5T.  ■. 

+  SJ>  +  ‘5JW  +  ^“‘al  +  a;’ +  '  3V  ’ V)  +  *  Vy) 

Again,  standard  fluid  dynamic  notations  are  used,  including  the  axial  velocity  u, 


the  radial  velocity  v,  the  temperature  T,  the  pressure  p,  the  thermal  conductivity 
Jt,  the  molecular  viscosity  p,  and  the  total  velocity  vector  V .  The  divergence  in 
cylindrical  coordinate  is  defined  by 


V  •  V 


Id.  du 

y  ay  ^  c 


The  total  energy  t  in  two  dimensions  is 


9  -  P£+  \p[u7  +  v7)- 

To  close  the  problem,  the  perfect  gas  relation  p  =  pRgT  is  also  required. 
When  written  in  vector  notation,  Eq.  (3.1)  becomes 

dQ  dE  dF  dEv  dFv 

dt  dz  dy  dx  dy 

in  which  dependent  variables  are  included  in  the  vector  Q  defined  by 

Q  =  y(p,p*,pv,e)T , 


(3.2) 


(3.3) 
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convection  terms  are  expressed  by  the  inviscid  flux  vectors  E  and  F  defined  by, 


’  pu 

pv 

E  =  y 

pu 3  +p 

,  F  =  y 

puv 

n 

puv 

pv  +  P 

.(«  +  p)u. 

.  (e  +  p)v . 

and  the  vectors  Ev  and  Fv  contain  the  second  order  viscous  diffusion  terms, 


Ev  =  y 


Mife-ift) 

<*(fc  +  ft) 

U(ft  +  ft)-M«(tft-!ft)  +  *ftJ 


(3.5) 


r  0 

F  -  ''•ft  + 

F°~v 

.««(!;  +  s?) + i*  - 1 1*) 

The  source  vector  H  includes  all  source  terms  associated 
etry  and  all  remaining  viscous  terms, 


(3.6) 


with  axisymmetric  geom- 


H~  +  ■  (3,7) 
-  “I£(MUv)  -  §^(MV2)  . 

In  this  form,  the  corresponding  equations  in  planar  two  dimensions  can  be  easily 
obtained  by  dropping  H  and  setting  the  j/’s  in  Eq.  (3.3)  through  Eq.  (3.6)  to  unity. 

To  facilitate  numerical  computation  over  arbitrary  geometries,  Eq.  (3.2)  is 
transformed  to  a  general  coordinate  system  by  a  transformation  defined  by 


t  -  £{x,y) 


n  =  n(*,v) 

where,  £  and  rj  are  usually  aligned  with  the  streamwise  and  the  cross-stream  direc¬ 
tion,  respectively. 


44 


The  transformed  equation  takes  the  form, 


dQ  dE  dF  -  dEv  dFv 
dt  '  dx  +  dy  *  dx  +  dy 


(3.8) 


in  which  the  strong  conservative  form  is  achieved  by  placing  all  the  metric  terms 
(£*»£v,etc)  inaide  the  derivatives  and  cancelling  the  arising  terms  by  the  metric 
identity.  The  new  dependent  variable  Q  is 

Q  ~  j(p,pu,pt/,e)T 

where,  J  is  the  Jacobian  of  the  coordinate  transformation,  which  is  defined  by, 


Inviscid  flux  vectors  in  the  general  coordinate  system  become 


pU 

pV  • 

puU  +  ftp 
PVU  +  ZyP 

£  y 

' 

puV  +  r)Tp 
pvV  +  rjvp 

L  (e  +  p)U  . 

i  (<  +  p)V  J 

in  which,  U  and  V  represent  contravariant  velocities  in  the  general  coordinate  sys¬ 
tem, 


U  =  +  Zvv 

V  as  T)x\l  +  T)VV. 

Viscous  flux  vectors  can  be  expressed  as 

fr  _  £zEv  +  £y  Fv 

Ev-  j 

r  _  VtF\i  +  HyFu 

K-  J 

while  the  source  vector  is  simply  fl  =  H/J. 

For  typical  high  Reynolds  number  flows,  a  highly  stretched  grid  is  required  in 
the  normal  direction  near  the  wall  in  order  to  resolve  the  large  normal  gradient  inside 
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the  boundary  layer.  Aa  a  consequence,  the  grid  size  in  the  atreamwise  direction 
generally  cfcnnot  be  refined  enough  to  resolve  corresponding  streamwise  gradients 


due  to  the  limitation  on  computer  storage  even  if  streamwise  diffusion  terms  are 
retained  in  the  complete  Navier-Stokes  equations.  Therefore,  for  high  and  even 
fairly  low  Reynolds  number  flows,  streamwise  diffusion  terms  can  manytimes  be 
neglected  without  losing  accuracy.  The  resulting  equation  set  is  referred  to  as  the 
thin-layer  Navier-Stokes  (TLNS)  equations  [llj, 


dQ  dE  dF  -  dFu 
dt  +  di  +  dr,  ~  +  dr, 


(3.9) 


The  TLNS  equations  will  be  used  as  the  governing  equations  for  viscous  calculations 


in  this  study. 

After  the  thin-layer  approximation,  the  viscous  flux  vector  Fv  and  the  source 
vector  H  become 

0 

+  asf^ 

ajf*  +  «3 

0 

p- 

-3  -  §»7| _ 

where, 

«i  =  +  til) 

M 

os  =  -r,tr,  v 

03  =  n[nl  + 

nrfc ,  3  2\ 

a<  =  +  nv) 

and  Cp  is  the  constant  pressure  specific  heal.  _'he  viscous  term  dFv/dr]  can  be 


further  rearranged  a? 


dK 

d’l 


d_ 
dr ] 


(*• 


dQx 

dr. 


+  R  2 


dQj 

dr] 


(3.10) 
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,/ 


in  which,  R i  and  i?j  are  4  x  4  matrices  defined  by 


1 

'0 

0 

0 

O' 

r  0 

0 

0 

0  * 

Ri  = 

0 

0 

a3 

Q3 

03 

0 

0 

*3  = 

0  0 

0 

0 

0 

0 

0 

0 

.0 

0 

0 

0. 

.04 

3 

an  —  a 

3 

03. 

The  vectors  Q  j  and  Qi  are  defined  by, 

»  tf*  , 

=  «,*/,«),  Qi  -  (e/p,u5,vs,uu)  . 


In  this  form,  the  viscous  dissipation  in  the  energy  equation  is  separated  from  the 
remaining  viscous  terms  so  that  formulations  for  which  viscous  dissipation  can  be 
neglected  are  easily  obtained  by  setting  £3  to  zero.  A  further  advantage  of  this 
splitting  is  that  matrices  R j  and  Ri  contain  only  metric  terms  and  properties  of 
the  gas  (viscosity,  thermal  conductivity  and  specific  heats).  For  cases  where,  n,k 
and  Cp  are  constants  or  nearly  so,  this  division  makes  the  linearization  of  the  viscous 
terms  particularly  easy.  For  turbulent  flows  when  these  quantities  vary  rapidly,  this 
form  separates  them  from  the  dependent  variables,  making  it  possible  to  identify 
their  effects  on  convergence  more  clearly. 

With  the  substitution  of  Eq.  (3.10)  into  Eq.  (3.9),  we  readily  have 


£2  +  ££  ,  ££  _  a  .  L(r  ££i  ,  „  *g», 

dt  di  dr]  dq  1  di]  2  dt]  ' 


(3.11) 


This  form  will  be  used  for  discussion  in  the  remaining  part  of  this  chapter.  Note 
that  the  Euler  equations  for  axisymmetric  two-dimensional  inviscid  flows  can  be 
obtained  by  dropping  all  the  viscous  diffusion  terms  in  Eq.  (3.11),  this  results  in 


dQ  dE  dF  - 

IT  +  a?  +  a?  =  "• 

where,  the  vector  Hi  represents  the  inviscid  source  vector, 


(3.12) 
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The  TLNS  approach  described  above  differs  from  the  classical  boundary  layer 
theory,  in  which  both  streamwise  diffusion  and  normal  pressure  gradient  are  ne¬ 
glected  as  a  result  of  the  order  of  magnitude  analysis.  Three  major  advantages  of 
the  TLNS  approach  over  the  classical  boundary  layer  theory  are: 

1.  Flows  containing  large  normal  gradients,  such  as  thick  boundary  layer  or  small 
reverse  flow  regions  can  be  readily  computed  by  using  the  TLNS  equations, 
since  the  pressure  gradient  is  retained  inside  the  boundary  layer. 

2.  The  pressure  inside  tne  boundary  layer  couples  with  the  pressure  variation  in 
the  inviscid  core  region  in  the  TLNS  approach,  hence  the  pressure  distribution 
along  the  cross-stream  direction  can  be  sot  ed  in  a  coupled  fashion  without 
need  for  the  inviscid-viscous  patching,  which  is  traditionally  used  in  the  classical 
approach  for  the  numerical  solution  of  high  Reynolds  number  nozzle  flows  {4,5j. 

3.  For  internal  flows,  the  TLNS  equations  automatically  conserve  mass,  while  the 
inviscid-viscous  calculation  based  on  the  boundary-layer  approach  generally  ig¬ 
nores  the  mass  inside  the  boundary  layer.  Although,  the  effect  of  this  mass 
deficit  on  the  inviscid  flow  is  compensated  by  offsetting  the  wall  contour  ac¬ 
cording  to  the  displacement  thickness,  the  patching  procedure  does  not  contain 
proper  mass  flow  rate.  This  is  especially  critical  for  flows  with  thick  boundary 
layers. 

3.2  The  Implicit  ADI  Scheme 

As  indicated  in  Chapter  1,  the  well-developed  alternating  direction  implicit 
(ADI)  scheme  has  been  extensively  used  to  solve  compressible  gas  dynamics  prob¬ 
lems.  In  this  section,  this  technique  is  formulated  for  the  axisymmetric  two- 
dimensional  TLNS  equations. 
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The  Euler  implicit  algorithm  for  Eq.  (3.11)  can  be  formally  expressed  as 


n+t 


=  0. 


(3.13) 


A  local  truncated  Taylor  series  expansion  can  be  used  to  linearize  all  the  flux  vectors 
and  the  source  term.  For  example,  E  is  linearized  by  using 


En+l  =  En  +  AAQ 


(3.14) 


where  AQ  =  Qn+I  -  Qn  and  A  is  the  Jacobian  matrix  dE/dQ  given  by 

0  Zv  0 

ts$-uU  t/  —  (nr  —  2)^«u  zvu  -  (-y  -  l)(*v  (n-l)(* 

ZyQ-vU  Uv  -  (ri  -  i)Zvu  u-('i-2)tvv  (n-i)^ 

U( 2*-n*)  UhLp  -  *)  -  h  -  l)uf/  ^w(nr^  —  4»)  —  (nr  —  l)vC/  iV 

in  which,  $  is  defined  by  4>  =  ^^-'(u2  +  v2). 

Similarly,  vectors  F  and  H  can  be  linearized  as 

/’n+1  =  Fn  +  BAQ  (3.15) 

Hn+l  =  Hn  +  DAQ  (3.16) 


where  B  and  D  are  Jacobian  matrices  dF/i)Q  and  dH jdQ.  The  matrix  B  can  be 
obtained  by  replacing  V  with  V  and  £  with  r?  in  the  A  matrix  given  above.  The 
Jacobian  matrix  D  is 

■  0  0  0  0  - 

p  _  djj  0  ^23  0 

<^31  d  32  ^33  d  34 

■di i  2  d< 3  0  . 
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where, 


2  t]x  d  uv  J 

dl'  =  '  y' 

3  J  dr)Kp  y 

#  4  uv  2  t)T  d  u  J 

U31  =  —  +  —  - — -  —  -u - (—  •  — ) 

y  3  py2  3  J  dr]''  p  y 

A  I  ,xU  2  fjr  8  1  A 

<*32  =  -(if  -  1)-  +  -M—  r-(-  ■  -) 

y  3  J  drt  p  y 

v  Alt 

<*33  =  -b-l)----^ 
y  3  py2 


<*34  = 


1  ~  1 


a  -  4  r/r  <3  4  A 

4‘  3J3r,(  p  'y1  3;^(  p  '  yJ 

2  7i  t  d  pv  J 

42  3  J  dr\  p  y 

43  3  J  dr)[  P  '  y’  3  J  dr,  p  y' 

The  viscous  terms  can  be  linearized  according  to 

n+  1 


dQ  i  dQi  d  d 

(3.17) 


where,  only  Q\  and  Q 2  are  linearized.  The  viscous  Jacobian  matrices  are  defined 
by  Bu\  =  OQ\/dQ  and  Bv2  =  dQi/dQ,  they  are  found  to  be 


Bvi  =  - 


1 


u  1 

p  p 


0  0  0 
0  0 


0 


0 


0  0  1 


and 


By,  2  =  — 


e 

7 


o  o 


i  -i 

p 

0  0 


~2~  2- 

p  (■ 

~2<Jj  0  2K  0 

I  -2—  *  ^  0 

L  p  n  n 


Substituting  Eq.  (3.14)  through  Eq.  (3.17)  into  Eq.  (3.13),  we  have 
3  d  d  3  3 

{I  -  MD  +  At(  —  A  -r  —  B  -  —{Ri~By,i  +  R2  —  Bu2)\}AQ  - -&tR  (3.181 
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where,  the  residual  vector  R  is, 


9E  dF  k  d  <Rd(>'  +RdQ'\\ 

d^{Rl~3T  +  R^)] 


(3.19) 


All  spatial  derivatives  in  Eq.  (3.18)  imply  that  central-differencing  will  be  used  for 
the  discretization. 

The  solution  of  the  Euler  implicit  scheme,  Eq.  (3.18),  requires  the  solution  of 
a  high  band-width  block  matrix  for  each  iteration,  unlike  the  one-dimensional  case, 
where  the  matrix  is  only  block  tri-diagonal.  The  computer  storage  and  the  CPU 
time  required  in  solving  this  high  band-width  matrix  are  very  large  for  typical  two- 
dimensional  problems.  Hence,  a  procedure  like  Eq.  (3.18)  is  impractical.  The  ADI 
or  Approximate  Factorization  (AF)  algorithm  arises  under  this  situation.  The  basic 
idea  of  a  ADI  scheme  is  to  split  the  left  hand  side  operator  of  Eq.  (3.18)  into  two 
parts  as, 


{I-AtD  +  At-~A){I  -  AtD)'1 

\I~AtD  4-  At  —  B  —  At~{R\  — +  B2"Bv2))AQ  -  ~ &tR 


(3.20) 


The  first  operator  on  the  left  hand  side  of  Eq.  (3.20)  contains  only  £  derivatives,  and 
the  last  one  contains  only  r?  derivatives.  The  solution  of  Eq.  (3.20)  can  be  broken 
into  two  steps.  For  the  £  direction  sweep,  the  equation 


(/  -  AtD  +  Af-rrA)  AQ’  -  -A tR 

is  solved.  After  obtaining  A Q'  at  each  grid  point,  ',he  jj  direction  sweep  equation, 
d  d  d  d 

(/  -  AtD  +  A t  —  B  -  At  —  (R\  ~£vi  +  R2  —  Bv2)\AQ  =  (/  -  AtD)AQ‘ 
orj  dr\  dt]  arj 

is  solved  for  A Q  over  the  entire  domain.  The  dependent  variable  is  then  updated 
according  to 

Qn+i  =  Qn  +  AQ. 
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Each  sweep  now  requires  only  the  solution  of  a  block  tri-diagonal  matrix.  This 
approach  saves  CPU  time  and  computer  storage. 

3.2.1  Boundary  Condition*  for  the  ADI  Scheme 


The  boundary  procedures  discussed  in  Section  2.2.1  can  be  easily  extended 
to  the  present  two  dimensional  calculation.  First  we  notice  that  the  similarity 
transformations 

A  =  M(A(Mf 1 ,  B  =  M„  A„M~ 1 

which  convert  Jacobian  matrices  into  diagonal  matrices  and  Ar,  exist  for  both 
Jacobian  matrices  A  and  B.  The  left  and  right  eigenmatrices  for  the  Jacobian 
matrix  A  are 


Mi  = 


m: 1  = 


u 

pk2 

V 

-pki 

L^T  p{uk2~vkx) 

1-0 

c* 

h-  1 

-Aciu  +  k,  v 

P 

Ai 

p 

(9,0 

0,0 

— •  and  k2  — 

"7fc  + 

iv  n 

1 

jTc  Tfc 


(7-1);*  -V 


^ _ pc 

~  v/2(7-1)c+  v/2(7-  1) 


9  —  k\u  +  k2v. 

The  eigenmatrices  Mn  and  A/~  1  can  be  obtained  by  the  substitution  of  £  with  r; 
in  Mi  and  M ^ 1 .  respectively.  The  transformed  matrices  A <  and  A,,  are  diagonal 
matrices  given  by 

\(  =  diag(U,  U,U  +  C^V  -  C(),  A„  =  diag(V,V,V  +  C„,V'  -  Cn) 
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in  which,  Ct  =  yj £\  +  £\c  and  =  yjnl  +  »/*c.  The  non-vanishing  elements  of 
A(  and  A„  are  the  eigenvalues  of  A  and  B,  respectively. 


For  demonstration  purposes,  the  boundary  procedures  for  the  flow  through  a 
nozzle  are  discussed.  Extension  to  other  types  of  flows  is  straightforward. 


As  indicated  earlier,  streamwise  diffusion  is  negligible  for  a  high  Reynolds  num¬ 
ber  flow.  Hence,  inviscid  MOC  boundary  conditions  are  applicable  at  the  inlet  and 
the  exit  plane.  If  the  inflow  is  supersonic,  all  eigenvalues  are  positive,  thus  four 
boundary  conditions  have  to  be  specified  in  order  to  determine  four  unknowns  at 
the  inlet.  In  other  words,  the  dependent  variables  must  be  completely  specified.  If 
the  inflow  is  subsonic,  only  the  fourth  eigenvalue  is  negative  since  U  <  C(.  This 
implies  three  specified  boundary  conditions  together  with  one  decoupled  character¬ 
istic  equation  must  be  imposed  at  the  inlet.  One  reasonable  choice  is  to  the  specify 
stagnation  pressure  P°,  the  stagnation  temperature  T°,  and  the  flow  angle  v/u. 
If  the  selection  matrix  L~  is  chosen  a a  L~  =  diag[0,0,0, 1)  and  the  vector  fl  is 
defined  as  0  ^  (P°,  T°,  v/u,  0),  a  boundary  procedure  similar  to  Eq.  (2.14)  can  be 
obtained  by  multiplying  Eq.  (3.20)  by  L~M{  and  combining  the  resulting  equation 
with  dU/dQ. 


No  special  treatment  is  needed  at  the  exit  if  the  outflow  is  supersonic.  For 
subsonic  outflow,  the  back  pressure  is  usually  specified.  Thus,  a  procedure  similar 
to  Eq.  (2.15)  can  be  applied  at  the  exit. 

At  the  centerline  of  the  nozzle,  the  dependent  variable  Q  vanishes,  since  y  is 
identically  zero.  Therefore  AQ  is  always  zero  at  the  centerline.  Flow  variables  can 
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be  calculated  after  each  iteration,  by  symmetrical  conditions, 


v  =  0 

r=° 

dy 

*£=0 

dy 

dT  n 

A  ~ 

dy 


At  the  wail,  the  MOC  type  procedure  can  be  obtained  for  inviscid  flows  by 
specifying  V  =  0  to  reflect  the  tangency  condition  and  selecting  the  decoupled  char¬ 
acteristic  equations  corresponding  to  outgoing  characteristics.  For  viscous  flows,  the 
MOC  procedure  is  not  valid,  the  no-slip  conditions, 


u  =  0 
v  =  0 


(3.21) 


together  with  zero  normal  pressure  gradient  and  the  specified  wall  temperature. 

(3.22) 


*£=0 

dn 


T  -  Tu 


can  be  used  instead.  Here,  Tw  is  the  specified  wall  temperature  and  n  denotes  the 
direction  normal  to  the  wall.  The  last  equation  can  be  replaced  by  adiabatic  wall 
or  specified  heat  flux  conditions  with  ndT/dn  =  q"{t). 

There  are  two  different  methods  to  apply  these  no-slip  conditions  at  the  wall. 
First,  &Q  can  be  set  equal  to  zero  at  the  *vaii  (13]  when  solving  the  discretized 
equation,  Eq.  (3.20).  Then  flow  variables  are  calculated  according  to  Eqs.  (3.21) 
and  (3.22)  by  using  the  updated  Q  from  interior  nodes.  In  this  way,  solutions  at  the 
wall  lag  those  of  the  interior  nodes  by  one  iteration  step,  thus  we  can  refer  to  this 
method  as  an  explicit-type  boundary  procedure.  Alternatively,  the  unknowns  at  the 
wall  can  be  coupled  to  the  unknowns  at  ♦he  interior  points  |14]  by  relating  &Q  at 
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wall  to  those  at  interior  nodes  according  to  the  no-slip  conditions.  In  this  approach, 
the  flow  variables  are  forced  to  obey  Eqs.  (3.21)  and  (3.22)  at  the  new  time  level. 
This  approach  solves  the  unknowns  at  wall  and  interior  nodes  simultaneously,  and  is 
referred  to  as  an  implicit-type  boundary  procedure.  In  the  following,  both  boundary 
procedures  are  implemented  at  wall  boundaries  and  their  effects  are  identified. 
3.2.2  Stability-Analysis  of  the  ADI  Scheme 

The  stability  behavior  of  the  implicit  ADI  scheme  has  so  far  not  been  iden¬ 
tified.  By  applying  the  double  Fourier  transform[24)  to  the  fully  implicit  scheme, 
Eq.  (3.18),  we  have 

I|Qn+1  =  LiQn 

so  that  the  amplification  matrix  G  is 

G  =  LJlL2 


where,  L\  -  I  +  Cri  and  Lj  —  I  with  Cpi  given  by 


_  .  At  At  At 

Cpi  =  -  A tD  +  i — -  sinw/A  +  i- — sinw„B  +  2 - r(l  -  cosu/n)(/?i  Bv\  +  j??u2 ) - 

A£  At7  An 

For  the  implicit  ADI  scheme,  L\  and  Lj  are  L\  ~  I  +  C>/  +  Cadi  and  Li  - 
l  +  Cadi ,  with 


Cadi  -  - 


Af^ 


sin u £  smu> r,A(I  -  .\tD)~lB 


+  2- 


A  t‘ 


sinu^(cosw„  -  1  )A(I  -  Af£>)"‘  {R\ Bvi  ~  R2Bv2). 


A^  An' 


The  matrix  Cadi  is  the  contribution  due  to  the  approximate  factorization.  In 
these  expressions,  uie  and  u/n  represent  the  wave  numbers  in  £  and  n  directions, 
respectively. 

As  stated  earlier,  the  numerical  stability  is  controlled  by  the  magnitude  of  the 
maximum  eigenvalue  of  the  matrix  G.  The  variation  of  the  maximum  eigenvalue 
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as  and  change  from  0  to  n  constitutes  a  three-dimensional  surface.  To  make 
the  results' more  readable,  only  the  variation  along  the  diagonal  line  =  u^)  on 
the  spectral  plane  is  calculated.  Figures  13  and  14  compare  the  results  of  the  algo¬ 
rithm  described  above  with  and  without  approximate  factorization.  The  flow  Mach 
number  is  0.5  and  the  Reynolds  number  is  104  for  both  cases.  The  time  step  size 
At  is  determined  by  the  CFL  number  based  on  U  +  C(.  The  maximum  eigenvalue 
is  shown  to  be  always  less  than  unity  for  the  fully  implicit  scheme,  and  hence  is 
unconditionally  stable  for  the  two-dimensional  TLNS  equations.  As  the  CFL  num¬ 
ber  increases,  the  maximum  eigenvalue  at  moderate  wave  numbers  (around  n/2) 
decreases.  This  implies  that  the  convergence  rate  will  be  speeded  up  by  increasing 
the  CFL  value. 

For  the  ADI  scheme,  the  approximate  factorization  alters  the  shape  of  stabil¬ 
ity  curves  substantially,  especially  when  CFL  is  large.  In  general,  the  maximum 
eigenvalue  at  the  mid-wavenumber  condition  is  much  higher  than  that  of  the  fully 
implicit  scheme  and  approaches  unity  as  the  wavenumber  increase.  The  results  on 
the  figure  suggest  that  the  ADI  algorithm  has  an  optimum  convergence  rate  at  a 
finite  value  of  CFL. 

Effects  of  the  flow  Mach  number  are  indicated  in  Fig.  15,  where  consecutive 
cases  are  calculated  for  several  Mach  numbers  with  a  CFL  number  of  10  for  the  ADI 
scheme.  The  eigenvalues  approach  unity  for  a  very  low  Mach  numbei  ( 1 0“ 2 ) ,  thus 
slow  convergence  can  be  expected.  As  the  Mach  number  increases,  the  eigenvalues 
near  low  and  high  wave  numbers  decrease,  thus  better  convergence  is  implied.  Al¬ 
though  not  shown  here,  effects  of  the  flow  Reynolds  number  on  stability  curves  a'e 
less  prominent  as  compared  to  those  of  the  flow  Mach  number  and  the  CFL  number 
for  the  present  local  stability  analysis. 
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igure  13.  Stability  results  of  2-D  fully  implicit  scheme  without 
approximate  factorization 


Maximum 


Figure  15.  Stability  results  of  2-D  implicit  ADI  scheme  with 

approximate  factorization  for  various  Mach  numbers 


The  Bow  through  an  axisymmetric  converging-diverging  nozzle  with  the  wall 
shape  defined  by  the  radius  variation, 

r(x)  -  -2{AR  -  l)x3  +  3 ARx1  +  1 

is  calculated  for  illustration.  Here  AR  is  the  ratio  of  the  throat  radius  to  the 
inlet  radius.  The  geometry  of  the  nozzle  and  a  representative  grid  are  shown  in 
Fig.  16.  Computational  results  are  shown  in  Figs.  17-19.  Two  typical  cases  are 
calculated,  including  a  pure  subsonic  Bow,  and  a  transonic  How  through  the  c-d 
nozzle  with  AR  =  0.8.  The  grid  is  50  x  30  with  30  in  the  cross-stream  direction  for 
inviscid  calculations  and  is  50  x  50  with  a  strong  clustering  near  the  wall  for  viscous 
calculations.  The  viscous  grid  is  shown  in  Fig.  16. 

The  computational  results  have  indicated  that  the  resulting  Mach  numbers 
at  the  entrance  are  around  0.2  and  0.4  for  the  subsonic  and  the  transonic  cases, 
respectively.  Figure  17  shows  the  L-2  norm  of  the  change  in  the  dependent  variable, 
AQ/Q,  versus  the  number  of  iterations  for  the  inviscid  calculations  by  using  the 
inviscid  grid.  Corresponding  results  for  the  viscous  calculations  based  on  the  viscous 
grid  are  shown  in  Fig.  18.  Both  inviscid  and  viscous  results  show  that  the  transonic 
case  converges  faster  than  the  pure  subsonic  case.  This  is  due  to  the  low  Mach 
number  effects  of  the  subsonic  case. 

The  comparison  of  Fig.  17  and  Fig,  18  indicates  that  ;nviscid  calculations  con¬ 
verge  faster  than  corresponding  viscous  calculations.  Ther*  are  two  possible  reasons 
for  this;  one  is  the  viscous  diffusion  inside  the  boundary  layer,  the  other  is  tne  grid 
stretching  near  the  wall.  To  understand  which  was  controlling,  the  inviscid  calcula¬ 
tions  are  done  on  the  viscous  grid  for  both  subsonic  and  ,ransonic  cases.  The  results 
are  also  shown  on  Fig.  18  ^s  is  seen,  inviscid  calculations  based  on  the  viscous 
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Figure  16.  Nozzle  geometry  for  2-D  transonic  calculations 
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grid  give  almost  the  same  convergence  rates  as  viscous  calculations.  Therefore,  we 
can  conclude  that  the  grid  stretching  deteriorates  the  convergence. 

The  optimum  CFL  number  for  all  calculations  is  around  5,  which  is  typical 
for  a  factorization  scheme  At  a  CFL  number  as  low  as  the  order  of  10,  further 
numerical  experiments  with  the  ADI  scheme  show  that  the  explicit  or  the  implicit 
wall  boundary  procedure  has  no  effect  on  the  convergence.  Computed  flowfield 
results  are  shown  in  Fig.  19  where  Mach  number  contours  for  the  viscous  transonic 
case  are  plotted. 

The  ADI  algorithm  has  also  been  applied  to  a  pure  supersonic  calculation.  To 
avoid  the  presence  of  shock  waves,  a  15  degree  conical  nozzle  with  an  expansion 
ratio  (the  exit  area  to  the  inlet  area)  of  30  and  an  inlet  Mach  number  of  1.02 
is  calculated.  Figure  20  shows  convergence  curves  for  both  viscous  and  inviscid 
results.  Comparisons  with  Fig.  17  and  Fig.  18  show  that  the  inviscid  supersonic  case 
converges  faster  than  the  inviscid  transonic  calculation,  while  the  viscous  supersonic 
case  is  as  slow  as  the  viscous  subsonic  case.  Although  these  preliminary  supersonic 
calculations  show  fairly  good  results,  further  numerical  experiments  on  a  higher 
expansion  ratio  nozzle  (100:1)  have  indicated  that  the  analysis  code  based  on  the 
central-differenced  ADI  scheme  is  difficult  to  start  with  an  arbitrary  assigned  initial 
guess.  Usually  a  very  low  CFL  number  (around  1)  has  to  be  used  at  the  initial 
stage.  Also,  the  convergence  rate  deteriorates  as  the  area  ratio  increases.  Additional 
experiments  of  transonic  flows  through  the  convergent-divergent  nozzle  described 
above  show  that  as  the  size  of  the  divergent  section  increases  by  a  certain  amount, 
the  central-differenced  code  becomes  unstable  and  diverges  and  the  upwind-central 
differencing  method  described  later  is  recommended. 

The  computational  results  above  illustrate  that  the  implicit  ADI  procedure  is 
efficient  for  flows  at  transonic  speeds.  For  flows  which  are  predominantly  super- 
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Figure  20,  Convergence  of  2-D  ADI  scheme,  inviscid  k  viscous  calculations, 
supersonic  flows 
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sonic,  such  as  those  flows  through  high  expansion  ratio  nozzles,  the  implicit  ADI 
procedure  is  inefficient  and  sometimes  even  unstable.  By  contrast,  our  stability  re¬ 
sults  have  indicated  that  high  Mach  rumber  flows  should  converge  faster  than  low 
Mach  number  flows.  The  reason  for  this  contradiction  is  as  follows.  Although  high 
Mach  number  flows  have  stronger  damping  effects  than  lower  Mach  number  flows, 
once  the  flow  becomes  supersonic,  the  equation  set  is  hyperbolic  tn  the  streamwise 
direction  since  all  eigenvalues  of  the  Jacobian  matrix  A  are  positive.  This  hyper- 
bolicity  implies  only  upstream  events  can  affect  downstream  events,  however,  the 
central-difference  formulation  in  the  AD!  scheme  allows  upstream  propagation.  In 
the  earlier  stages  of  a  computation,  this  upstream  influence  keeps  on  propagating 
wrong  information  (from  the  unconverged  solution)  from  downstream  back  to  up¬ 
stream,  and  consequently  slow.*  down  the  convergence.  This  also  explains  why  the 
coda  is  difficult  to  start  for  supersonic  cases. 

3.3  Stability  Consideration  of  Unwind  Algorithms 

As  we  have  seen  in  the  previous  section,  the  imp'icit  ADI  algorithm  becomes 
inefF.v  ient  when  the  flow  is  predominantly  supersonic  due  to  tho  fact  ‘.bat  the  central- 
difference  formulation  allows  the  downstream  influence  of  the  unconverged  rolulioi. 
10  propagate  upstream.  To  avoid  this  unwanted  upstream  propagation,  upwind 
schemes  appear  to  be  desirable.  The  basic  idea  of  a  upwind  difference  i.'  »o  model 
th<*  physics  correctly  by  using  a  difference  stencil  which  involves  only  windward 
information,  hence  “wrong"  propagation  is  prohibited.  The  upwind-difference  for¬ 
mula!  io  n  is  chosen  to  develop  suitable  numerical  algorithms  for  supersonic  calcula¬ 
tions.  Based  on  the  successful  application  in  one-diinersional  flows,  che  flux-vetior 
.ipi.Liing  algorithm  ( 16}  will  again  be  utilized  to  develop  unwind  algorithms,  be¬ 
fore  formulating  the  detail  of  the  numerical  procedures,  the  stability  analysis  of  a 
modeled  scalar  equation  is  considered  for  several  possible  solution  procedures 
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To  simulate  present  two-dimensional  viscous  flows,  the  Burger’s  equation, 


du  .du  du  du  d2u 

It  +a  Tx*“  Tx  +  %  =  "a? 


(3.23) 


is  chosen  as  the  modeled  equation,  in  which  a*  and  a~  are  positive  and  negative 
constants,  respectively.  The  second  and  the  third  terms  on  the  left  hand  side  are 
used  to  simulate  the  flux-splitting  in  the  streamwise  direction,  where  the  flow  is 
predominantly  supersonic.  Again,  these  terms  must  be  differenced  according  to  the 
windward  directions.  Although  second  order  upwind  differencing  will  be  used  to 
formulate  the  numerical  procedures  for  Navier-Stokes  calculations;  for  the  modeled 
equation,  only  first  order  accurate  differencing  is  considered.  The  last  term  on  the 
left  hand  side  irr. plies  that  central-differencing  will  be  retained  since  no  preferable 
windward  direction  exists  on  the  cross  plane.  The  second  order  term  on  the  right 
hand  side  is  used  to  model  the  thin-layer  viscous  diffusion. 

3.3.1  The  Fully  Implicit  Scheme 


Direct  application  of  the  Euler  implicit  scheme  to  Eq.  (3.23)  results  in 

d  d  d  ^2 

(1  +  Af(a+  —  +  a-  —  +  6- - /j—v)]Au  -  -Atr 

dz  dz  dy  dy2 


(3.24) 


where,  r  is  the  residual  given  by 

( 


,  +  du  _du  ,du  d2u  ' 

r  =  (a  Tx+a  d-x+bTy-“W] 


The  amplification  factor  gn  for  Eq.  (3.24)  can  be  easily  found  to  be 


1 


:fi  = 


l  +  C>/ 


where,  Cri  is 


Cjr;  =  C;(  1  -  cos uz  + » sin u>x)  +  i}~  (cos  ut  -t- »  sin ut  -  1 )-»-  t<7v  sin  +  2^(1  -  cos  vv). 
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The  parameters,  <r+,<r~,and  oy  are  CFL  numbers  defined  by 


a~  At 
Ax  ’ 


bAt 
Ay ' 


The  parameter  u  is  the  von  Neumann  number  defined  by 


v  = 


ix  At 

Ay7 


As  stated  earlier,  the  direct  solution  of  a  fully  implicit  scheme  is  impractical 
in  multi-dimensions,  hence  approximate  factorization  is  necessary.  A  number  of 
approximate  procedures  can  be  identified  for  solving  the  left  hand  side  operator  in 
Eq.  (3.24).  Of  them,  three  factorization  procedures  will  be  discussed. 


3.3.2  The  Standard  API  Scheme 


The  first  approximate  procedure  is  to  split  the  operator  in  Eq.  (3.24)  as 


[i  +  At(a  +  ~~  +«'—)!(!  +  At(b^  -  »-*~)\Au  =  -  Air. 


(3.25) 


Equation  (3.25)  is  analogous  to  the  implicit  ADI  procedure,  except  upwiuding  has 
been  used  in  the  streamwise  direction.  This  splitting  generates  error  terms  on  the 
left  hand  side  of  order  At7.  The  resulting  amplification  factor  qAd\  for  Eq.  (3.25) 
is 


9  ADI  = 


1  +  ^ADl 

1  +  Cri  +  Cadi 


where, 


C ADI  —  "  cos  w,  -r  i  sin  uz)  sin  uv  4-  \a ~ ov{cz»u;x  +  t  sin  uiz  -  1 )  sin  u\ 

-f  2a * i/(l  -  cosw,  +  » sin wt)(l  -  cost* >y) 

4-  2o~  i/(cos  uiz  4-  t  sin  ur  -  1)(1  -  cot  u>„) 
is  the  summation  of  factorization  error  terms. 
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3.3.3  The  Diagonally  Dominant  ADI  Scheme 

i 

The  second  splitting  under  consideration  is  the  diagonally  dominant  alternating 
direction  implicit  (DDADI)  method  suggested  by  Lombard  (14).  This  splitting, 
when  first  order  upwind  differencing  is  used,  can  be  expressed  as 


.  Ata+  A1/,  d  d2  Ata 

|d-— u.-1,,+AI(»— |  d+ 


d  92 

Al  = -A,r 


(3.26) 


where  d  is  the  summation  of  diagonal  elements  given  by  d  =  1  +  ~  °-)- 

Equation  (3.26)  can  be  solved  directly  by  using  alternating  forward  and  backward 
sweeps.  Alternatively,  it  can  also  be  solved  by  using  the  line  Gauss-Seidel  relaxation 
by  MacCormack  [15]  and  Chakravarthy  ( 17] .  This  method  has  also  been  referred 
to  as  the  single-level  scheme  by  Lombard  (19],  and  the  LU  scheme  by  Yoon  and 
Jameson  jl8j. 

The  line  relaxation  method  for  the  DDADI  splitting  is  discussed  as  follows, 
First,  Eq.  (3.26)  can  be  split  into, 


I* + =  -A,r + fr* j 


(3.27) 


and 


.  ,  .  „  d  d2  . .  .  . .  .  At  _  A 

\d  +  At(6- - /i—  Au  s  dAu  -  —a  Aul+i,; 

ay  ay 2  Ax 


(3.28) 


Equation  (3.27)  can  be  rearranged  as 


|d+4*(fcA_MjL)]Au*  =  -A<|a+ 


u  —  u  , 
Ax 


+  (a 


du  du 
—  + 
dx  dj 


d 1  v 


(3.29) 


Equation  (3.29)  is  lower  bi-diagonal  and  can  be  solved  by  marching  from  the  up¬ 
stream  toward  the  downstream  at  each  cross-stream  station  after  obtaining  u'_  , 
from  the  iteration  at  previous  stations.  Similar  to  this  forward  marchi-g  proce¬ 
dure,  we  can  define  a  symmetrica  backward  marching  procedure  by  substituting 
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Au*  from  Eq.  (3.27)  into  Eq.  (3.28)  and  neglecting  unnecessary  terms  in  order  in 
make  the  resulting  right  hand  side  residual  to  be  consistent  with  the  steady  state 
equation,  r  =  0,  as  Au  is  driven  to  zero.  The  resulting  backward  marching  equation 
can  be  written  as, 


-  U"+V,y  u«.j 

Ax 


+  du  du  fl5u  ‘ 
(a + 


(3.30) 


Equation  (3.30)  now  can  be  solved  by  marching  from  the  downstream  toward  the 
upsirsam  using  the  updated  value  of  uj1^1 1  ■  at  the  t  +  1  station.  The  combination  of 
Eq.  (3.79)  and  Eq.  (3.30)  completes  one  iteration  step  of  the  line-relaxation  method 


for  the  DDADI  splitting. 


The  amplification  factor  g *  for  the  intermediate  forward  marching,  Eq.  (3.29), 


is 


.  _  _ 1  -  ^“(cosu*!  +  isincj!) _ 

u"  1  +<j/(1  -  cosw,  +  isin-j*)  -  aZ  +  ioyn\t\uy  +  2u(l  -  coswy) 

and  the  amplification  factor  g "  for  Eq.  (3.30)  is 

••  _  un+>  _ _ 1  +  (cos u>x  -  t  sin uz ) 

u*  1  i-  at  ■+•  (cos u,  +  s  sinwz  -  1)  +  iov  ainu/v  -f  2u{l  -  co3u/v) 

The  overall  amplification  factor  is  then 

9DDAD '  =  (3.31) 


3.3.4  The  Parabolized  ADI  Scheme 

The  third  algorithm  considered  is  the  splitting, 

d  $  ^  2  ^ 

11  +  At(a*ib+  bTy  +  =  "Af.  (3-3*) 

As  we  shall  see  in  the  next  section,  the  first  operator  on  the  LHS  is  similar  to  a 
parabolized  Navier-Stoke  (PNS)  marching  operator,  hence.  Eq.  (3.32)  is  referred  to 
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as  a  P NS- ADI  splitting.  The  amplification  factor  for  this  splitting  can  be  found  to 


be, 


9PNS-ADI 


1  “I-  CpNS 
l  +  Cri  +  Cpus 


(3.33) 


in  which, 

Cpns  =  1  —  co8wj  +  *  sinwx)(cosu;x  +  tsinu/x  —  1)  + 

ia~Oy(co9u>t  +  isinu/x  —  l)sinwy  +  2o~i/(cosu>x  +  »sinu;x  -  1)(1  -  cosu>y). 

Note  that  if  a~  vanishes,  the  error  Cpss  is  zero  and  the  PNS-ADI  splitting  is 
equivalent  to  the  fully  implicit  scheme. 


3.3.5  Algorithm  Comparisons 

The  amplification  factors  for  all  four  algorithms  noted  above  can  be  numerically 
computed  over  the  entire  spectral  plane  for  wx  and  u>y  ranging  from  0  to  ir.  Again, 
to  make  results  more  concise  and  readable,  only  the  variation  along  the  diagonal 
line  (u/z  =  u >y)  on  the  spectral  plane  is  calculated.  Two  cases  are  studied,  they  are, 

1.  Subsonic  :  o+  =  —o~  =  ov  =  u  =  CFL 

2.  upersonic:  o+  —  oy  =  v  —  CFL ,  a~  —  0 

The  first  case  is  analogous  to  the  subsonic  flow,  where  both  upstream  and 
downstream  propagations  are  significant.  The  second  case  simulates  supersonic 
flow  by  setting  a~  =  0.  Both  cases  assume  uniform  CFL  numbers  in  the  x  and  y 
directions. 

Results  of  these  two  cases  are  plotted  on  Fig.  21  and  Fig.  22  for  a  CFL  of 
10.  For  the  first  case  both  the  fully  implicit  and  the  DDAD1  algorithms  give 
monotonically  decreasing  eigenvalues,  and  these  eigenvalues  reach  their  minimums 
at  the  wavenumber  n.  The  eigenvalue  first  decreases,  then  increases  to  a  local 
maximum  near  the  wave  number  n  for  the  ADI  and  PNS-ADI  schemes,  In  the 
second  case,  all  algorithms  except  the  ADI  splitting  give  monotonically  decreasing 
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eigenvalues.  This  it  due  to  the  factorization  errors  of  the  ADI  splitting  even  when 
a~  is  zero.  ' The  eigenvalue  of  the  DDADI  splitting  is  slightly  smaller  than  that  of 
the  fully  implicit  scheme  and  the  PNS-ADI  splitting  in  the  second  case.  Again,  the 
naturally  dissipative  characteristic  of  upwind  schemes  is  observed  by  noting  that 
the  eigenvalues  are  much  less  than  unity  near  the  wave  number  it. 

As  a  final  test  case,  the  effects  of  CFL  numbers  on  the  magnitude  of  the  eigen¬ 
values  are  shown  in  Fig.  23  for  the  subsonic  case  ( a~  ^  0).  As  can  be  seen,  increasing 
the  CFL  number  tends  to  magnify  the  eigenvalues  for  the  PNS-ADI  splitting,  thus 
a  CFL  number  of  103  is  inferior  to  a  CFL  number  of  10  as  far  as  convergence  is 
concerned.  On  the  other  hand,  the  DDADI  splitting  is  insensitive  to  the  CFL  num¬ 
ber.  Effects  of  the  CFL  number  on  ADI  and  fully  implicit  schemes  are  similar  to 
those  on  PNS-ADI  and  DDADI  schemes,  respectively. 

Among  the  four  algorithms,  the  fully  implicit  scheme  and  the  DDADI  splitting 
can  be  expected  to  give  better  convergence  rates  than  the  other  two  splittings, 
since  they  have  smaller  eigenvalues  over  the  entire  spectrum.  The  DDADI  splitting 
is  superior  to  the  other  two  approximate  procedures  and  is  as  good  as  the  fully 
implicit  scheme.  The  PNS-ADI  splitting  will  give  better  convergence  than  the 
standard  ADI  splitting  since  they  have  similar  behavior  in  subsonic  regions  (a~  *  0) 
and  the  PNS-ADI  scheme  is  more  dissipative  in  supersonic  regions  (a"  =  0).  All 
four  algorithms  can  be  shown  to  have  eigenvalues  always  less  than  unity,  hence  they 
are  unconditionally  stable  for  the  two-dimensional  Burger's  equation. 

Based  on  the  results  of  stability  analysis,  the  proposed  hybrid  algorithm  which 
uses  upwind  difference  in  the  streamwise  direction  and  central  difference  in  the 
cross-stream  direction  proves  to  be  unconditionally  stable  for  the  Burger’s  equation 
for  fully  implicit  and  all  three  approximate  procedures  investigated. 


Eigenvalue 
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The  stability  results  discussed  above  are  based  on  first-order  upwind  differ- 
encing  and  the  scalar  modeled  equation.  For  algorithms  based  on  more  accurate 
second-order  upwind  differencing  and  the  vector  equation,  the  stability  results  may 
be  quite  different  from  the  results  given  above.  For  instance,  the  approximate  fac¬ 
torization  error  term  associated  with  the  product  of  a+djd x  and  a~d/dx  in  the 
Parabolized  ADI  scheme  (Eq.  (3.32))  reduces  to  zero  for  the  vector  equation  since 
A*  A~  is  identically  zero  for  both  the  subsonic  and  supersonic  cases.  This  implies 
that  the  PNS-ADI  procedure  will  give  better  convergence  than  that  predicted  by 
the  stability  analysis  based  on  the  Burger’s  equation.  In  the  following  discussion, 
this  hybrid  algorithm  will  be  utilized  to  formulate  numerical  algorithms  for  viscous 
supersonic  calculations. 

3.4  Algorithms  for  Viscous  Supersonic  Flows 

As  is  well  known,  the  hyperbolic  nature  of  supersonic  flows  allows  inviscid  prob¬ 
lems  to  be  computed  in  a  single  pass.  This  capability  is  lost  in  viscous  flows  for  two 
reasons.  The  first  is  the  streamwise  diffusion,  but  this  is  weak  for  high  and  even 
moderate  Reynolds  number  flows  and  may  frequently  be  neglected,  which  results  in 
the  TLNS  equations.  The  second  and  more  significant  reason  is  the  subsonic  region 
near  the  wall,  where  upstream  influences  exist.  From  a  physical  point  of  view,  this 
subsonic  region  allows  information  to  propagate  upstream.  This  upstream  propa¬ 
gation  prohibits  a  single  pass  solution  and  renders  the  space-marching  procedure 
unconditionally  unstable.  In  the  PNS  approach,  this  upstream  influence  is  removed 
by  a  parabolized  procedure  which  neglects  the  contribution  associated  with  this 
upstream  influence  so  that  a  single  sweep  solution  is  allowed  (26-29].  For  problems 
with  significant  upstream  influence,  such  as  thick  boundary  layers  inside  the  nozzles 
1 5 1  and  separated  flows,  global  iterations  are  required,  The  main  focus  in  this  chap¬ 
ter  falls  into  this  category,  thus  upstream  effects  must  be  retained  by  considering 
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the  global  iteration  from  the  Navier-Stokes  equations.  The  simplification  to  PNS 
procedures  will  be  discussed  in  the  next  chapter. 

Based  on  the  results  of  stability  consideration,  in  discretizing  the  governing 
equations,  we  can  take  advantage  of  the  predominantly  supersonic  character  of  the 
flowfteld  by  using  flux-vector  splitting  in  the  streamwise  direction,  while  retain¬ 
ing  central  differencing  in  the  cross-stream  directions.  Again,  using  Euler  implicit 
differencing  in  time,  the  result  can  be  expressed  as 


Q’+'-Q”  ,SE*  ^  £  SCl^^ 

At  d£  d£  dr)  dr)  1  dr)  2  dr) 


n+l 

)]  =0(3.34). 


If  we  assume  the  flux  vector  E  is  homogeneous,  it  can  be  split  according  to 


£  =  £++r  =  A+Q  +  A-Q. 


Matrices  A+  and  A~  are  the  split  Jacobian  matrices  given  by 

A+  =  M(A+M{  \  A~  =  A/<A'A/"' 

in  which, 

At  =  (A<  +  |A<|)/2,  A"  =  (A*  -  |A<|)/2. 

The  quantities  E+  and  E~  have  to  be  upwind  differenced  according  to  the  signs 
of  their  individual  eigenvalues.  In  supersonic  regions,  E~  vanishes  and  Eq.  (3.34) 
provides  an  algorithm  which  can  be  solved  by  a  marching  procedure  in  the  £  direc¬ 
tion.  The  difficulty  is  that  the  supersonic  region  is  coupled  to  the  subsonic  region 
where  upstream  influence  exists.  Consequently,  both  regions  must  be  solved  multi¬ 
ple  times.  The  following  discussion  formulates  numerical  procedures  for  Eq.  (3.34) 
based  on  the  algorithms  considered  in  Section  3.3. 
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3.4.1  The  Standard  ADI  Procedure 


The  first  procedure  considered  is  the  standard  approximate  factorization  |8,9) 
of  the  upwind  differenced  system.  In  this  procedure,  the  £  and  rj  directions  are 
split,  and  the  discretized  version  of  Eq.  (3.34)  in  delta  form  becomes, 


\l-AtD  +  +  ~AT)\(1  -  AID)'  ' 

I/-AID+  Al^-B  -  +  R,~ B  ,)|A(j  =  -AtR' 

1  drj  dr)  dr)  dr) 


where,  the  residual  vector  PJ  is, 


R'  = 


dE'  dF 
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dQj 
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(3.35) 


(3.3G) 


The  Jacobian  matrices  A?  and  Af  are, 
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dQ  ' 
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dE~ 

dQ 


where  the  subscripts  t  are  used  to  distinguish  the  true  Jacobi&ns  Af  from  the 
matrices  A± .  The  first  operator  in  Eq.  (3.35)  is  block  penta-diagonal  for  second 
order  upwind  differencing  and  is  block  tri-diagonal  for  the  first  order  scheme.  The 
solution  procedure  of  Eq.  (3.35)  is  similar  to  that  o'  fhe  implicit.  ADI  scheme  and 
will  not  be  repeated  here. 


3.4.2  The  Diagonally  Dominant  ADI  Procedure 


The  second  approximate  factorization  procedure  is  the  DDADI  method.  To 
express  this  algorithm,  we  must  discretize  the  equation  in  £  and  t  before  factoring. 
The  explicit  discretization  in  p  is  not  necessary  to  specify  the  algorithm  and  the 
derivatives  in  rj  imply  central  differencing.  The  basic  philosophy  of  the  DDADI 
procedure  is  to  place  as  many  terms  as  possible  on  the  diagonal  element  before 


splitting.  This  procedure  can  be  expressed  as 


{/  -  A<i>  -  A,-  +  All 


|/-A(D  +  *|i(^  -yt,-)]"' 
.  At 


</  -  a«j> + *£**♦  +  A<i^r  +  £*  -  £(*,£*.,  +  ac 

=  -At/?' 

(3.37) 

in  which,  /c  =  1  4-  *c/2.  The  quantity  k  is  zero  for  first  order  upwind  differencing 
and  is  one  for  second  order  upwind  differencing. 

Similar  to  the  derivations  discussed  in  Section  3.3.3,  tue  DDADI  splitting  will 
be  solved  by  the  line  Gauss-Seidel  relaxation  method,  which  includes  the  forward 
marching, 

<B'+A<||-B  -  A(fl,  i-j},,  +  = 


dr] 


..as:.,)"- 

al  A£  +  2A* 

and  the  symmetrical  backward  marching, 

<£>'+A![|j£>  -  ±{R,±B„  +  B,  Aflrt)|)A(j  = 


(3.38) 


-  Af{ 
.dE+ 


(^u)n+1-fer  _ 

Af  2A£ 


n+  1 


-  5  ,  dQ{  dQi  x 
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(3.39) 

where,  A(j*  =  Q‘  —  Qn  and  AQn+1  =  Qn+I  —  Q'  are  used  to  update  variables 
immediately  after  each  forward  or  backward  sweep,  respectively.  The  diagonal  term 
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Both  forward  and  backward  marching  operators  involve  the  solution  of  a  block 
tri-diagonal  matrix.  The  combination  of  Eq.  (3.38)  and  Eq.  (3.39)  complete*  one 
iteration  step. 


3.4.3  The  PNS-ADI  Procedure 


The  third  approximate  factorization  procedure  considered  is  also  a  forward- 
backward  scheme  that  is  designed  to  give  a  PNS-like  sweep  in  the  forward  direction 
and  a  partial  backward  sweep  that  is  required  only  in  those  regions  where  the  flow 
is  subsonic.  Because  of  its  analogy  with  PNS  algorithms,  this  procedure  is  referred 
to  as  a  PNS-ADI  scheme.  This  scheme  can  be  described  as, 


{/  -  At D  +  At\~-Af  +  ~D-  --(/?,  +  *2|-0u2)j)(/  -  AtD)~' 

oq  oq  oq  on 


(I  -  AtD  +  At  —  A;)AQ  =  -AtR' 


(3  40) 


As  will  be  shown  later,  the  first  operator  is  exactly  the  time-iterative  version  of  a 
PNS  procedure,  as  wili  be  discussed  in  detail  in  chapter  4.  Again,  this  operator 
is  block  tri-diagonal  at  each  (  location.  The  last  operator  reduces  to  the  identity 
operator  in  jupersonic  regions  where  vanishes,  and  is  only  block  bi-diagonal 
in  subsonic  regions.  Thus,  in  supersonic  flows,  the  PNS-ADI  procedure  becomes  a 
marching  procedure,  but  in  subsonic  regions  it  retains  the  influence  of  the  upstream 
acoustic  waves.  Because  the  backward  operator  is  only  necessary  inside  the  subsonic 
layer,  the  computational  time  involved  in  one  iteration  for  this  PNS-ADI  algorithm 
is  less  thar  that  for  the  DDADI  algorithm. 


3.4.4  Direct  Solution  by  Flowfleld  Partitioning 


For  a  typical  high  Reynolds  number  supersonic  flow,  the  flowfleld  can  almost  be 
solved  by  a  pure  marching  procedure.  Only  the  thin  subsonic  layer  adjacent  to  the 
wall  prevents  this  marching  algorithm.  This  suggests  a  direct  solver  based  on  the 
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flowfield  partitioning  according  to  the  distinct  physical  nature  of  this  predominantly 
supersonic  flow.  As  is  given  in  Fig.  24,  the  flowfield  is  partitioned  into  three  parts. 
Domain  A  is  composed  of  all  the  points  that  are  supersonic.  Domain  B  includes  all 
the  points  that  are  subsonic  (although  it  may  also  include  some  supersonic  points). 
Domain  C  is  the  interface  line  between  domain  A  and  domain  B.  It  is  chosen  such 
that  all  points  on  the  line  (domain  C)  are  supersonic.  Let  QAt  Qb ,  and  Qc  denote 
all  the  dependent  variables  in  these  domains.  These  unknowns  are  related  by 

AaQa  +  CaQc  —  Ra 

BbQb  +  CbQc  —  Rb  (3.41) 

Ac  Qa  +  BcQb  +  CcQc  =  Rc 

which  indicates  that  the  unknowns  in  domain  A  are  coupled  only  to  unknowns  at 
domain  C  and  are  independent  of  unknowns  at  domain  B.  Similar  statement  holds 
for  domain  B,  while  domain  C  is  coupled  to  both  A  and  B. 

Since  all  points  are  supersonic  at  domain  A,  the  operator  Aa  can  be  efficiently 
solved  by  a  marching  procedure.  The  operator  Bb  is  block  tri-diagonal  with  block 
size  equal  to  4  x  Jb,  (Jb  denotes  the  number  of  points  in  r?  direction  inside  domain 
B).  This  block  tri-diagonal  matrix  is  large,  but  the  size  decreases  rapidly  as  the 
number  of  subsonic  points  decreases  (as  for  thin  subsonic  layers).  Both  QA  and 
Qb  can  be  solved  easily  if  the  unknowns  Qc  are  given.  In  fact,  Qc  can  be  sohed 
according  to 

\AcAa'Ca  4-  BcBj,1  Cp  -  Cc\Qc  =  AcAa  1  IiA  +  Be  Bp  '  Ru  ~  Re  (3  12) 

which  is  obtained  by  the  substitution  of  the  first  and  the  second  equations  in 
Eq.  (3.41)  into  the  lost  equation.  The  left  hand  side  operator  of  Eq.  (3.42)  is  a 
dense  matrix  with  the  size  equal  to  4  x  /,  with  /  representing  the  number  of  points 
in  the  (  direction  for  the  entire  computational  domain.  This  matrix  can  be  effi¬ 
ciently  solved  by  an  iterative  method  in  4  to  6  iterations,  because  the  influence  of 
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downstream  points  on  upstream  points  is  weak  if  all  points  along  domain  C  are 
supersonic.'  This  partitioning  technique  allows  a  direct  solution  procedure  to  be 
obtained.  The  direct  solution  without  any  approximate  factorization  implies  the 
CFL  number  can  be  as  large  as  possible,  thus  rapid  convergence  can  be  obtained. 
However,  the  computational  time  involved  in  one  iteration  for  this  direct  method  is 
much  more  than  those  for  approximate  procedures  if  too  many  subsonic  points  are 
involved. 

3.5  Algorithm  Comparisons 

Based  on  the  physics  of  supersonic  flows,  four  different  upwind  algorithms  asso¬ 
ciated  with  their  stability  analyses  are  discussed  in  the  last  two  sections.  However, 
the  computational  efficiency  of  these  numerical  procedures  is  yet  to  be  identified. 
Before  comparisons  to  be  made,  two  concepts  must  be  paid  attention  to.  First, 
we  note  that  all  four  methods  described  above  are  concerned  with  using  different 
algorithms  to  solve  the  same  equations.  The  residuals  on  the  right  hand  side  of 
Eq.  (3.35),  Eq.  (3.37),  Eq.  (3.40),  and  Eq.  (3.41)  are  identical  in  terms  of  both  par¬ 
tial  differential  equations  and  finite  difference  representations.  As  A Q  is  driven  to 
zero,  all  four  methods  provide  the  same  steady  state  solution  as  given  by  Eq  (3.36) 
A  check  of  the  converged  solutions  from  the  computer  codes  verifies  that  these  so¬ 
lutions  are  identical  to  six  or  seven  digits.  Second,  since  all  methods  give  the  same 
solutions,  the  only  thing  to  compare  is  the  path  a  specific  algorithm  takes  to  steady 
stale.  Therefore,  it  is  of  interest  to  compare  both  the  number  of  iterations  required 
for  A Q  to  reach  the  machine  accuracy  and  the  total  CPU  time  required.  The  num¬ 
ber  of  iterations  required  shows  the  numerical  efficiency  of  each  algorithm,  while  the 
CPU  time  required  indicates  the  cost  per  iteration  step.  The  latter  is  afTected  by 
the  arithmetic  operations  involved  in  an  algorithm  and  the  computer  architecture. 
Present  results  are  for  the  scalar  machine,  IBM3090-180.  Slightly  different  results 
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may  be  obtained  for  a  vector  machine. 

The  test  problem  for  the  comparison  follows  is  the  laminar  supersonic  flow 
through  a  15°  conical  nozzle  with  an  expansion  ratio  of  30.  The  nozzle  geometry 
and  the  21  (axial)  by  40  (radial)  grid  are  shown  on  Fig.  25.  At  the  starting  plane, 
a  constant  Mach  number  flow  (M  =  1.02)  with  zero  contravariant  velocity  V  was 
chosen.  The  nozzle  Reynolds  number  based  on  the  throat  radius  is  10s.  The 
converged  solution  showed  that  with  this  highly  stretched  grid,  the  subsonic  portion 
of  the  boundary  layer  had  grown  to  nine  points  at  the  exit  plane.  All  calculations 
are  done  without  artificial  viscosity  for  the  upwind-central  differencing  formulations. 
Boundary  conditions  are  implemented  by  the  procedures  given  in  Section  3.2.1. 

Convergence  rates  baaed  on  the  number  of  iterations  for  the  four  algorithms 
mentioned  above  are  shown  in  Fig.  26.  The  linear  convergence  on  these  semi- 
logarithmic  plots  until  machine  accuracy  is  reached  gives  indication  that  the  codes 
are  error-free.  All  calculations  started  with  initial  conditions  which  were  obtained 
from  the  single  forward  sweep  through  the  flowfield  with  the  PNS  algorithm  (The 
details  of  PNS  algorithms  will  be  given  in  chapter  4).  The  result  shown  for  each 
algorithm  corresponds  to  the  optimum  CFL  for  this  scheme.  These  optimums  are 
shown  in  the  figure.  Boundary  conditions  are  implemented  by  the  implicit  wail 
boundary  procedure.  In  terms  of  number  of  iterations  required  the  ADI  scheme 
is  seen  to  be  the  slowest.  It  also  has  the  lowest  optimum  CFL  of  5.  The  DDADI 
scheme  is  the  fastest  of  the  three  approximate  methods  and  converges  almost  as 
rapidly  as  the  direct  method.  The  DDADI  algorithm  converged  most  rapidly  with 
CFL  at  5000,  above  this  value,  convergence  becomes  independent  of  CFL. 

Figure  26  also  shows  that  the  PNS-ADI  algorithm  gives  excellent  convergence 
(it  converges  to  machine  accuracy  in  70  iterations).  This  rapid  convergence  was 
not  properly  predicted  by  the  stability  analysis  based  on  the  scalar  equation.  The 
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optimum  CFL  .'or  ‘he  PNS-ADI  procedure  is  around  100,  beyond  that,  convergence 
rates  tend  to  alow  down.  The  optimum  CFL  for  the  direct  solver  was  about  10° 
and  the  convergence  remained  the  same  for  CFL  numbers  up  to  1010.  The  conver¬ 
gence  rate  for  the  direct  procedure  was  somewhat  disappointing  in  that  it  required 
25  iterations  to  reach  the  machine  accuracy.  With  the  fully  implicit  scheme,  the 
numerical  procedure  can  normally  converge  to  machine  accuracy  in  8-  11  iterations 
(42|.  This  normal  convergence  was  not  obtained  with  the  current  analysis  code. 

Comparisons  of  convergence  rates  in  terms  of  CPU  time  are  shown  in  Fig.  27. 
The  ADI  method  is  seen  to  require  the  most  time  to  reach  machine  accuracy.  The 
direct  method  is  just  slightly  faster  than  the  ADI  method;  however,  if  we  could 
attain  a  factor  of  three  improvement  expected  above,  it  would  be  competitive  with 
the  fastest  procedure.  The  PNS-ADI  and  the  DDADI  schemes  are  quite  competitive 
and  are  about  a  factor  of  ten  faster  than  the  standard  ADI  scheme.  Although,  in 
terms  of  number  of  iterations,  the  PNS-ADI  method  is  slower  than  the  DDADI 
method,  the  deficit  is  seen  to  be  offset  by  less  computational  work  involved  in  one 
iteration  for  the  PNS-ADI  procedure. 

As  demonstrated  in  Chapter  2,  the  use  of  the  true  Jacobian  has  significant  ef¬ 
fects  on  convergence.  The  same  comparison  was  made  for  all  four  methods  by  using 
approximate  Jacobians  A+  instead  of  A?  on  the  left-hand  side  of  corresponding 
discretized  equations.  The  convergence  rates  in  terms  of  both  number  of  itera¬ 
tions  and  CPU  time  are  shown  in  Fig.  28  and  Fig.  29.  All  procedures  except  the 
ADI  method  slow  down  substantially  due  to  the  approximation  in  Jacobians.  The 
convergence  rate  (both  in  terms  of  iterations  and  CPU  time)  for  the  ADI  method 
actually  improves  slightly.  This  shows  that  the  convergence  rate  for  the  ADI  scheme 
is  insensitive  to  the  Jacobian  matrices.  Note  also  that  corresponding  optimum  CFL 
numbers  for  three  approximate  methods  change  if  approximate  Jacobians  are  used. 
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Figure  27.  Convergence  in  terms  of  CPU  time  of  2-D  supersonic  algorithms 
with  true  Jacobians 
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Figure  28.  Convergence  in  terms  of  number  of  iteration  of  2-D  supersonic 
algorithms  >  ith  approximate  Jacobians 
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Figure  29.  Convergence  in  terms  of  CPU  time  of  2-D  supersonic  algorith 
with  approximate  Jacobians 
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3.6  Nozzle  Flowfleld  Predictions 


To  demonstrate  the  capabilities  of  present  upwind-central  difference  algorithms, 
the  flowfields  in  two  nozzle  geometries  were  computed.  The  first  geometry  was  the 
15°  conical  nozzle  shown  on  Fig.  25  except  here  a  21  x  70  grid  was  used.  The  second 
nozzle  was  the  contoured  geometry  shown  in  Fig.  30,  with  a  grid  of  75  x  50.  The 
area  ratio  of  this  nozzle  was  272  :  1.  In  both  nozzles,  a  throat  radius  of  10  mm.  was 
chosen. 

The  properties  of  a  typical  rocket  nozzle  combustion  gas  were  used  for  all 
calculations,  including  -7  =  1.24  and  Cp  =  3043 J/ Kg°K.  Both  nozzles  were  run  at 
tne  stagnation  pressures  of  35  and  3.5  atm.,  corresponding  to  nozzle  throat  Reynolds 
numbers  of  1.4  x  10*  and  1.4  x  105.  The  stagnat'on  temperature  was  chosen  to  be 
3500 °F.  The  results  for  the  conical  nozzle  were  for  laminar  flows  with  the  molecular 
viscosity  varying  according  to  the  Sutherland  law, 


3/2  Tr  +  5 
T  +  S 


where,  ^ r  is  the  reference  viscosity  at  a  reference  temperature  of  TV  and  5  is  the 
Sutherland  constant.  Those  results  for  the  contoured  nozzle  were  calculated  for 
both  laminar  and  turbulent  flows.  For  turbulent  calculations,  the  algebraic  model 
by  Baldwin  and  Lomax  (12}  was  used. 

3.0.1  Verification  of  Solution  Accuracy 


To  verify  the  accuracy  of  present  upwind-central  difTerencing  algorithms,  the 
results  obtained  from  present  analysis  codes  are  compared  to  those  from  the  MOC 
procedure  (4).  Figure  31  plots  Mach  number  contours  for  the  inviscid  supersonic  flow 
through  the  high-expansion  contoured  nozzle.  The  upper  half  shows  results  from 
the  MOC  procedure  and  the  lower  half  presents  those  from  present  algorithm- 
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As  shown  here,  the  Mach  number  contours  computed  by  using  current  upwind- 
central  differencing  schemes  agree  very  well  with  the  MOC  predictions,  except  for 
the  discrepancies  near  the  centerline,  where  the  MOC  procedure  fails  to  resolve  the 
symmetry  conditions  due  to  the  presence  of  a  weak  oblique  shock. 

Figure  32  shows  the  comparison  between  current  algorithms  and  the  MOC  pro¬ 
cedure  by  plotting  the  computed  wall  pressure  distribution  along  the  axial  direction. 
Again,  the  inviscid  results  from  current  algorithms  show  excellent  agreements  with 
those  from  the  MOC  procedure.  The  corresponding  wall  pressure  distribution  for 
viscous  calculations  with  a  Reynolds  number  of  1.4  x  10s  is  also  shown  on  the  fig¬ 
ure.  The  viscous  wall  pressure  is  slightly  higher  than  that  of  inviscid  results.  This 
illustrates  that  viscous  flows  expand  less  due  to  the  presence  of  the  boundary  layer. 

3.6.2  Effects  of  Downstream  Boundary  Conditions 

All  numerical  algorithms  mentioned  above  require  a  downstream  boundary 
condition  at  the  subsonic  part  inside  the  boundary  layer,  since  the  eigenvalue  U  ~C{ 
is  negative  if  U  <  C (.  Previous  researchers  have  usually  implemented  this  boundary 
condition  by  extrapolating  from  Inside  the  computational  domain.  Although  this 
does  simplify  the  numerical  procedures  at  the  downstream  boundary,  it  violates 
the  physical  conditions,  especially  when  the  boundary  layer  is  thick.  Further,  an 
extrapolated  boundary  condition  does  not  allow  (lowfields  to  respond  to  downstream 
pressure  changes  as  they  must  in  physical  situations.  In  the  case  of  nozzle  flows,  a 
high  back  pressure  will  cause  the  boundary  layer  to  thicken  and  the  flow  lo  separate. 

One  case  of  interest  is  when  an  exhaust  nozzle  is  operated  in  an  altitude  facil¬ 
ity  where  the  ambient  pressure  in  the  facility  is  only  approximately  matched  to  the 
nozzle  expansion  characteristics.  This  mismatch  can  provide  significant  differences 
between  test  stand  performance  and  eventual  performance  in  space  because  of  the 
effect  of  the  ambient  pressure  on  the  nozzle  boundary  layer  at  the  exit.  An  analo- 
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Figure  31.  Comparison  of  Mach  number  contours  for  MOC  and 
upwind/centra!  difference  algorithm 
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Figure  32.  Comparison  of  wall  pressure  distribution  for  MOC  and 
upwind /central  difference  algorithm 
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gous  but  much  more  severe  result  of  this  mismatch  is  that  flow  separation  orrurs 
inside  the  thick  boundary  layer  of  an  attitude  nozzle  when  it  is  operated  at  sea-level 
conditions.  The  examples  shown  later  will  refer  only  to  the  case  of  small  back  pres¬ 
sure  mismatches  and  will  not  extend  to  the  massive  separations  that  are  observed 
at  sea  level  conditions.  Calculations  of  these  effects  require  that  the  downstream 
boundary  conditions  be  correctly  implemented. 

The  downstream  boundary  conditions  for  the  present  ir  .^ligation  are  based 
on  MOC  procedures.  At  the  early  stages  of  iteration,  a  certain  number  of  points 
at  the  exit  plane  will  become  subsonic  due  to  viscous  diffusion.  The  back  pressure 
is  then  specified  for  these  subsonic  points  on  the  exit  plane.  The  procedure  used 
is  similar  to  that  in  Eq.  (2.15).  As  the  iteration  proceeds,  reverse  flow  will  appear 
inside  the  subsonic  layer  if  the  back  pressure  is  high  enough.  For  these  reentry  flows, 
three  eigenvalues  (U ,  U ,  and  U -C(]  become  negative(assuming  the  reentry  velocity 
is  less  than  the  sonic  speed)  and  standard  inflow  boundary  conditions  are  applied. 
This  implies  three  conditions  must  be  specified  from  outside  the  domain  while  one 
characteristic  equation  must  be  used  for  information  coming  from  inside  the  domain. 
These  three  specified  quantities  are  chosen  as  the  stagnation  pressure,  the  stagnation 
temperature,  and  the  flow  angle,  corresponding  to  the  ambient  conditions  Since 
the  external  surroundings  are  interpreted  as  being  at  rest,  the  stagnation  pressure 
is  taken  as  the  back  pressure.  For  those  subsonic  points  with  positive  (outflow) 
streamwise  velocities,  only  the  back  pressure  is  specified  and  three  equations  are 
used  in  agreement  with  traditional  outflow  boundary  procedures.  The  decision  as 
to  whether  inflow  or  outflow  boundary  conditions  are  used  depends  on  the  signs  of 
eigenvalues  at  each  grid  point  as  determined  from  the  previous  iteration  The  above 
procedure  for  reentry  flows  is  reliable  for  modestly  sized  recirculation  regions  but 
eventually  breaks  down  for  large  recirculation  regions  because  oscillations  between 
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inflow  and  outflow  occur. 

Results  of  a  series  of  calculations  in  which  the  pressure  at  the  subsonic  part 
of  the  outflow  boundary  was  specified  are  shown  on  Figs.  33-38.  Figure  33  shows 
Mach  number  contours  for  the  conical  nozzle  at  the  lower  Reynolds  number.  The  top 
plot  shows  the  results  by  using  extrapolation  boundary  conditions.  This  boundary 
condition  resulted  in  a  back  pressure  to  stagnation  pressure  ratio  Pb/ P°  of  about 
2.5  x  10-3.  The  high  Mach  number  gradient  near  the  wall  gives  an  indication  of  the 
boundary  layer  thickness.  The  middle  plot  shows  the  effect  of  raising  the  back  pres¬ 
sure  to  Pb/P°  =  5  x  10-3.  At  this  back  pressure,  the  boundary  layer  is  thicker  and 
a  small  recirculation  zone  presents  near  the  exit.  As  the  back  pressure  is  increased 
further,  this  recirculation  region  continues  to  grow  and  to  propagate  upstream  as 
indicated  by  the  bottom  plot  where  Pb/P°  =  7  x  10" 3 .  Corresponding  calculations 
with  back  pressures  below  2.5  x  10-3  showed  that  the  boundary  layer  at  the  exit 
accelerated  and  became  thinner.  The  wall  temperatures  for  these  calculations  are 
3000 °K. 

The  dotted  lines  in  Fig.  33  represent  the  Mach  number  contours  from  one¬ 
dimensional  inviscid  calculations  corresponding  to  the  same  stagnation  conditions. 
As  is  seen,  the  presence  of  the  boundary  layer  in  the  two-dimensional  calculations 
has  resulted  in  less  expansion  and  thus  slower  speeds  near  the  exit. 

As  a  further  indication  of  the  character  of  these  flows,  the  axial  velocity  (u) 
profiles  at  the  exit  plane  for  various  back  pressure  levels  are  shown  on  Fig.  31. 
This  figure  shows  how  the  boundary  layer  grows  as  the  back  pressure  is  increased 
and  the  rate  of  thinning  of  the  boundary  layer  as  the  back  pressure  is  decreased. 
Furthermore,  the  width  of  the  recirculation  zone  at  the  exit  plane  is  clearly  shown. 

The  results  for  turbulent  flow  in  the  contoured  nozzle  at  the  higher  Reynolds 
number  (P°  =  35a<m)  are  shown  in  Fig.  35.  Again,  the  top  plot  shows  Mach 
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Figure  34.  Velocity  profiles  at  nozzle  exit  plane  for  conical  nozzle  using 
various  back  pressures 


0.5  x  10  \  ^ 


0.7  x  10  3 


Radius  (meter) 


5000 


Figure  38.  Velocity  profiles  at  exit  plane  for  laminar  flow  in  272:1 
contoured  nozzle 
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number  contours  for  the  extrapolated  boundary  condition  case.  The  middle  and 
the  bottom 'plots  show  the  results  of  higher  back  pressures.  The  extrapolation  case 
here  corresponds  to  the  back  pressure  ratio  Pb/P°  of  0.33  x  10  3.  The  middle  plot 
is  for  a  back  pressure  of  0.7  x  10-3  while  the  lower  one  is  for  l.l  x  10~3.  Similar 
thickening  and  separation  of  the  boundary  layer  is  observed. 

Because  of  the  interchange  of  momentum  inside  the  boundary  layer  for  tur¬ 
bulent  flows,  a  larger  increase  in  Pb!  P°  is  required  to  obtain  the  the  same  de¬ 
gree  of  separation,  as  is  verified  by  comparing  the  laminar  flow  results  shown  on 
Fig.  36.  For  this  laminar  calculation,  the  extrapolated  boundary  condition  cor¬ 
responds  to  a  lower  back  pressure  { Pb/P°  =  0.29  x  10~3)  than  to  the  turbulent 
case,  because  the  much  thinner  boundary  layer  allows  more  expansion  to  be  accom¬ 
plished  in  the  nozzle.  Also,  the  laminar  calculations  show  a  larger  separation  region 
at  Pb/P°  =  0.7  x  10-3  than  the  turbulent  boundary  layer  at  Pt,/P°  s  1,1  x  10~3. 
The  u  velocity  profile  at  the  exit  plane  for  laminar  and  turbulent  boundary  layers 
are  shown  on  Fig.  37  and  Fig.  38,  respectively.  The  turbulent  case  has  a  thicker 
boundary  layer  and  steeper  gradients  at  the  wall  than  the  laminar  calculation. 

3-6.3  Effects  of  Reynolds  Number  and  Wall  Temperature 

Although  calculations  with  two  different  Reynolds  numbers  have  been  shown  ii. 
the  last  section,  some  additional  results  of  changing  Reynolds  numbers  are  demon¬ 
strated  here  along  with  comparisons  of  the  effect  of  wall  temperature.  Figur?  39 
shows  Mach  number  contours  for  the  conical  nozzle  at  a  Reynolds  number  of  1.4  x  IQ4 
but  with  a  lower  wall  temperature  of  300 °K.  The  upper  plot  is  for  extrapolation 
boundary  conditions,  and  the  lower  one  is  for  a  back  pressure  of  7  x  10' 3 ,  analogous 
to  the  bottom  plot  of  Fig.  33.  These  results  show  that  colder  wall  temperatures 
give  a  much  thinner  boundary  layer.  As  a  consequence,  a  smaller  separation  region 
is  observed  as  compared  to  the  hot  wall  boundary  layer. 
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Figure  40.  Converged  Mach  number  contours  for  high  Reynolds  number 
Extrapolated  boundary  conditions 
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The  effects  of  Reynolds  numbers  on  the  conical  nozzle  results  can  be  seen  by 
comparing  the  Mach  number  contours  for  a  stagnation  pressure  of  35  atm  in  Fig.  *10 
with  those  for  the  lower  Reynolds  number  case  with  P°  =  3.5  atm  on  Fig.  33.  Again, 
a  thinner  boundary  layer  is  obtained  for  the  higher  Reynolds  number  case.  Only 
the  results  based  on  the  extrapolated  boundary  conditions  are  given  in  Fig.  40. 

As  a  final  comparison,  results  for  the  contoured  nozzle  are  shown  on  Fig.  41  for 
the  low  Reynolds  number  case  {3.5  atm)  at  two  different  v,all  temperatures,  3000P/C 
and  300 °K,  based  upon  turbulent  boundary  layer  assumptions.  As  for  the  laminar 
case,  the  colder  wall  temperature  results  in  a  much  thinner  boundary  layer. 

The  effects  of  Reynolds  number  can  be  seen  by  comparing  the  results  in  Fig.  4 1 
with  those  for  the  turbulent  case  in  Fig.  35.  Again,  a  lower  Reynolds  number  causes 
a  dramatic  increase  in  the  boundary  layer  thickness.  At  this  low  Reynolds  number, 
a  check  of  the  maximum  eddy  viscosity  in  the  boundary  layer  profile  reveals  that 
the  flow  is  no  longer  “turbulent"  although  the  turbulence  model  is  still  included. 

3.6.4  pLCqupM  wall  Qgpling  gnd  plows 

As  we  have  seen  in  previous  sections,  the  wall  temperature  has  a  significant 
effect  on  the  nozzle  flowfield.  While  specified  temperature  or  heat  flux  boundary 
conditions  are  frequently  imposed  at  wall  boundaries.  In  viscous  problems,  rocket 
nozzle  walls  are  in  general  regeneratively  cooled  by  propellant  flowing  inside  the 
wall.  This  poses  a  problem  when  neither  the  heat  flux  nor  the  wall  temperature  are 
known  a  priori  but  both  must  be  completed  as  part  of  the  nozzle  flowfield  solution. 
In  the  present  section,  we  develop  a  coupled  method  for  solving  the  wall  cooling 
flow  along  with  the  nozzle  flow. 

Figure  42  shows  the  schematic  of  a  nozzle  surrounded  with  cooling  tubes.  The 
coolant  is  assume  to  flow  from  the  exit  toward  the  inlet.  To  simplify  the  problem, 
the  following  assumptions  are  made: 


contours  for  272:1  nozzle 


at  lower  Reynolds 


no 


1.  The  outer  wall  of  the  cooling  tube  is  insulated. 

2.  The  wfcll  between  the  coolant  the  nozzle  flow  is  thin  and  of  high  conductivity 
so  that  the  temperature  of  the  cooling  liquid  can  be  taken  equal  to  that  of  the 
wal^Tu,). 

3.  The  cooling  liquid  has  a  constant  specific  heat  (Cj). 

Let  yw  and  n  denote  the  radius  and  the  inward  normal  direction  at  an  arbitrary 
wall  location  ({,q).  Referring  to  Fig.  42,  the  energy  balance  of  the  control  volume 
at  a  wall  location  gives 

mC,Tw  +  q"dA  =  mC,(Tw  - 

where  dA  represent  the  surface  area  wetted  by  the  cooling  liquid.  If  the  nozzle  wall 
is  completely  surrounded  by  cooling  tubes,  the  quantity  dA  can  be  expressed  in  the 
general  coordinate  system  as 

dA  =  ^j^y/vi+ri'ldt. 


By  the  Fourier's  law  and  geometrical  relations,  the  heat  flux  q"  is 


sfnl  +  n'i  at  dn 


in  which,  k  is  the  thermal  conductivity  of  the  gas.  Substituting  expressions  for  q" 
and  dA  into  the  energy  balance  equation,  we  readily  obtain 
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(3.43) 


where  y  is  the  normalized  wall  radius  (y  =  yw/yt)  and  B\  is  the  non-dimensional 
Biot  number  defined  by 
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The  non-dimensional  Biot  number  represents  the  ratio  of  wall  cooling  to  heat  con¬ 
duction  in  the  gas. 

The  new  thermal  boundary  condition,  Eq.  (3.43),  can  be  coupled  in  an  implicit 
manner  to  the  discretized  governing  equations  provided  that  the  temperature  at  the 
downstream  end  of  the  wall  (the  inlet  temperature  of  the  cooling  liquid)  is  given. 
The  numerical  procedure  is  based  on  the  DDADI  scheme.  The  wall  temperature  at 
the  exit  is  fixed  at  the  given  coolant  inlet  temperature.  In  the  discretized  equation 
of  Eq.  (3.43),  the  derivative  dTw/d (  is  backward  differenced  since  the  coolant  flows 
from  downstream  to  upstream.  The  derivative  dTw/drj  is  one-sided  differenced 
and  is  coupled  to  the  unknowns  of  interior  points.  The  implicit  treatment  of  this 
discretized  boundary  equation  is  similar  to  the  boundary  procedure  discussed  in 
Section  3.2.1. 

Typical  Mach  number  contours  for  supersonic  flows  through  the  high  expansion 
ratio  nozzle  given  in  Fig.  30  by  using  this  wall  cooling  boundary  condition  are  shown 
in  Fig.  43.  The  top  plot  is  for  adiabatic  wall  conditions,  the  middle  and  the  bottom 
plots  are  for  B i  =  103  and  Bi  =  10®,  respectively.  The  inlet  temperature  of 
the  coolant  is  500° K  for  the  last  two  cases.  The  Reynolds  number  is  1.4  <  104 
and  the  flow  is  assumed  laminar  for  all  cases.  Dramatic  changes  in  the  boundary 
layer  thickness  and  the  flowfield  near  the  exit  when  this  more  appiopriale  cooling 
condition  is  incorporated  can  be  observed.  The  results  shown  are  for  demonstration 
only,  for  practical  applications,  the  parameter  Bi  should  be  calculated  according  to 
the  true  wetted  area  of  cooling  tubes  and  real  properties  of  the  coolant. 

3.0.5  Nozzles  with  Subsonic  Inflow 

So  far,  the  examples  of  supersonic  nozzle  flowfield  predictions  we  have  seen  start 
from  an  arbitrarily  given  Mach  number  distribution  at  the  inlet.  For  real  nozzles, 
the  flow  enters  the  diverging  section  with  a  non-uniform  Mach  number  distribution. 
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To  take  account  of  this  effect,  the  computation  must  begin  with  the  subsonic  section 
of  the  noztfe.  As  we  have  shown,  the  transonic  flow  through  the  nozzle  throat  region 
can  be  efficiently  calculated  by  the  implicit  ADI  scheme.  For  flowfield  computations 
over  realistic  nozzles  which  contain  a  short  converging-diverging  section  and  a  very 
long  diverging  section,  the  following  procedure  is  recommended.  The  nozzle  can  be 
segmented  into  two  parts.  The  first  part  contains  the  entire  converging  section  and 
a  small  portion  of  the  diverging  section.  The  divergent  portion  is  chosen  sufficiently 
large  to  ensure  the  flow  at  the  flow  at  the  last  few  rows  of  grid  points  is  supersonic 
except  for  the  boundary  layer.  The  implicit  ADI  algorithm  can  be  applied  to  this 
transonic  portion  and  the  resulting  flowfield  near  the  exit  ran  be  used  as  the  input 
for  subsequent  supersonic  calculations  for  the  remaining  part  of  the  nozzle.  The 
previously  described  supersonic  algorithms  can  then  be  easily  applied  to  flowfield 
computations  in  the  diverging  section. 

Typical  laminar  results  of  the  computation  over  the  contoured  nozzle  starting 
from  subsonic  inflows  by  using  the  procedure  above  are  shown  in  Fig.  44.  The  grid 
is  300  x  50  with  300  in  the  axial  direction  and  the  Reynolds  number  is  1.4  x  1C4. 
Comparisons  with  previous  results  using  constant  Mach  number  at  the  inlet  (Fig.  41) 
show  that  the  two-dimensionality  near  the  throat  has  only  minor  effects  on  the 
flowfield  results  for  this  typical  example. 

3.6.0  Verification  of  Global  Conservation 

It  is  generally  agreed  that  for  flows  with  discontinuities,  the  strong  conservative 
form  of  the  equations  plays  an  important  role  in  global  conservation.  For  flowfields 
that  do  not  contain  discontinuities,  the  fully  conservative  form  is  sometimes  assumed 
to  be  less  important.  The  primary  application  of  present  numerical  algorithms  are 
for  rocket  nozzle  flowfield  predictions.  To  accurately  predict  the  flowfield  and  noz¬ 
zle  performance,  global  conservation  is  of  great  importance;  however,  "good"  nozzle 
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designs  will  generally  be  free  from  shocks  and  it  might  be  assumed  that  the  non¬ 
conservative  equations  are  adequate.  Because  of  the  importance  of  accurate  global 
mass  conservation  and  to  demonstrate  the  necessity  of  using  the  conservative  form, 
mass  conservation  was  monitored  in  all  analysis  codes  throughout  the  course  of 
this  study.  For  better  understanding  of  the  necessity  of  the  strong  conservative 
formulation,  a  few  computations  with  the  weak  conservative  form  of  the  governing 
equations  were  also  done  for  typical  high  expansion  nozzle  flows  without  disconti¬ 
nuities.  The  weak  conservative  formulation  is  identical  to  Eq.  (3.9)  except  that  the 
metric  coefficients  are  left  outside  the  derivatives, 
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The  global  mass  conservation  results  for  this  weak  conservative  and  the  strong 
conservative  formulations  are  compared  in  Fig.  45  for  the  high  expansion  ratio 
contoured  nozzle.  The  75  x  50  grid  shown  in  Fig.  30  are  used  for  both  cases.  As 
shown  on  this  figure,  the  strong  conservative  form  maintains  the  mass  flow  rate 
error  within  1%,  while  the  weak  conservative  form  gives  a  maximum  mass  error  of 
about  30%.  In  fact,  for  all  the  cases  computed  to  date  with  the  strong  conservative 
form,  including  calculations  for  a  nozzle  with  expansion  ratio  as  high  as  70 0,  the 
maximum  mass  flow  rate  errors  have  been  maintained  below  1%. 

The  reason  for  the  failure  in  weak  conservative  formulation  is  because  it  does 
not  conserve  the  mass  in  its  finite  difference  representation  while  the  fully  conser¬ 
vative  form  does.  Therefore,  the  results  shown  in  Fig.  45  is  to  be  expected.  The 
above  results  demonstrate  that  even  for  flows  without  discontinuities,  the  strong 
conservative  formulation  is  necessary  for  maintaining  global  conservation. 


quasi-conservative  scheme,  error=28% 
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CHAPTER  4 


THE  APPLICATION  OF  TIME-ITERATIVE  SCHEMES  TO 
THE  PARABOLIZED  NAVIER-STOKES  EQUATIONS 


As  indicated  in  Chapter  1,  the  parabolized  Navier-Stokes  equations  have  been 
extensively  used  as  an  alternative  to  the  Navier-Stokes  equations  for  the  solution 
of  compressible  as  well  as  incompressible  viscous  flows  due  to  their  computational 
efficiency.  To  assess  this  efficient  numerical  procedure  and  place  it  in  a  unified 
context  with  the  present  Navier-Stokes  procedures,  the  parabolized  technique  is 
also  addressed  here.  In  Chapter  3,  time-iterative  numerical  procedures  based  upon 
the  predominant  physics  of  the  flow  were  formulated  for  the  solution  of  thin-laver 
Navier-Stokes  equations.  Starting  from  these  time-iterative  Navier-Stokes  equa¬ 
tions,  it  is  shown  in  this  chapter  that  the  parabolized  equations  can  be  obtained  as 
a  subset  of  the  Navier-Stokes  equations  by  means  of  flux  splitting  procedures.  The 
traditional  pressure  gradient  splitting  devised  by  Vigneron  is  shown  to  be  a  specific 
type  of  flux  vector  splitting  while  the  method  of  characteristics  based  splitting  devel¬ 
oped  for  the  TLNS  equations  in  Chapter  3  is  also  shown  to  split  the  streamwise  flux 
vector  into  parts  with  positive  and  negative  eigenvalues  which  can  also  be  treated 
in  a  “parabolized”  sense.  This  new  PNS  procedure  is  formulated  by  neglecting  the 
flux  vector  with  negative  eigenvalues.  The  computational  results  obtained  by  using 
both  the  classical  and  the  new  PNS  procedure  are  compared  to  those  obtained  with 
the  TLNS  algorithms  to  verify  the  accuracy.  The  flux  splitting  interpretation  of  the 
PNS  procedure  allows  a  well-behaved  global  iterative  PNS  procedure  to  be  defined 
based  on  the  PNS-ADI  method  presented  in  Chapter  3. 
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The  thin-layer  Navier-Stokes  equations  as  given  in  Eq.  (3.11)  will  be  used  as  a 
starting  point  for  the  present  approach.  We  first  split  the  streamwise  flux  E  into 


two  parts, 


E  =  E+  +t 


where  the  eigenvalues  of  the  Jacobians  of  E+  and  E "  are  positive  and  negative, 
respectively.  With  the  substitution  of  Eq.  (4.1)  into  Eq.  (3.11),  we  have 

i  +  +  (4.2) 

dt  d£  di  drr  dr)  dt)  dr) 

Note  that  the  splitting  given  by  Eq.  (4.1)  is  used  conceptually  to  indicate  a  general 
expression  of  flux-vector  SDlitting,  which  might  represent  the  Steger  and  Warming 
splitting  defined  in  Chapter  3,  but  which  could  also  represent  any  other  splitting 
procedures.  As  examples,  we  will  in  addition  to  the  Steger- Warming  splitting  also 
discuss  a  splitting  based  upon  Vigneron’s  [26]  parabolization  method. 

Although  the  approach  is  equally  applicable  to  homogeneous  or  inhomogeneous 
flux  vectors,  for  simp  Jty,  we  have  restrict  to  the  homogeneous  case  where  E  =  AQ 
with  A  ~  dE/dQ.  Flux-splitting  of  the  homogeneous  vector  is  then  reduced  to 


splitting  the  matrix  A  as, 


A  =  A+  +  A~ . 


Here  the  eigenvalues  of  A+  are  positive  and  those  of  A  are  negative.  From  the 
homogeneity  of  the  vector  E,  we  have 

E+  =  A  +  Q ,  E~  --  A~Q 

which  obviously  satisfies  Eq.  (4.1). 

In  principle,  both  the  streamwise  and  the  cross-stream  fluxes  could  be  split, 
instead,  we  follow  the  traditional  PNS  approach,  for  which  central  differences  are 
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generally  used  in  the  cross-stream  direction  while  upwina  differences  are  used  in 
the  streamwiie  direction.  Thus,  just  as  we  did  for  the  TLNS  equations  in  Chapter 
3,  we  have  split  only  the  streamwise  flux  h. 

Using  Euler  implicit  differencing  in  time,  the  discretized  version  of  Eq  (4.2) 
can  be  expressed  as 

d  d  d  d  d  d 

{I-AtD  +  Al[—At  +  —A-  +  ~B--f(Rl^-Bvl+R2^-Bv2)\)AQ  -  -A tR' 

di  di  '  dr)  dr)  dr)  dr) 

(4.4) 

where  R'  has  been  previously  defined  by  Eq.  (3.36).  Again,  ♦he  true  Jacobians  of 
E+  and  E~  are  indicated  by  A*  and  Af ,  respectively  and  the  spaxial  derivatives  in 
Eq.  (4.4)  must  be  treated  consistently  on  both  the  left  side  and  the  right  hand  side. 
Note  aiso  that  the  derivatives  containing  A*  and  Af  must  be  upwind  differenced 
in  the  manner  defined  previously.  Efficient  solution  of  Eq.  (4.4)  requires  cme  sort 
of  approximate  factorization  of  the  type  discussed  in  Chapter  3. 

41.1  Spotting  Based  on  Characteristics 

As  the  first  of  our  two  specific  flux  splitting  procedures,  we  begin  with  one 
baaed  on  the  method  of  characteristics.  As  we  have  seen  in  the  previous  discussion, 
the  matrix  A  can  be  diagonalized  according  to  the  similarity  transformation  given 
by 

A  =  Sr~'AM( 

where  the  subscript  of  A  has  been  dropped  for  simplicity.  The  matrix  A/e  is  com¬ 
posed  of  the  right  eigenvectors  of  the  matrix  A  end  has  been  previously  defined. 
The  diagonal  matrix  A  contains  four  entries,  U ,  U ,  U  +  C( ,  and  U  -  Ce.  A  straight¬ 
forward  splitting  suggested  by  Sieger  and  Warming  [10]  is, 

.  (A  +  lAll 

A  = - 

2 

(A-iAj) 


(1.5) 
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in  which,  |A|  refer*  to  the  matrix  composed  of  the  absolute  values  of  the  elements 
of  A.  FromEq.  (4.5)  we  readily  obtain 

=  A^A+M^1,  A~  =  (4.6) 

For  the  homogeneous  case,  the  split  flux  vectors  thus  become 

E+  =  A+Qt  E~  -  A~Q. 

According  to  Eq.  (4.5),  for  supersonic  flows,  the  matrix  A"  is  identically  zero,  and 
A+  is  equal  to  A.  For  subsonic  flows,  these  matrices  are  A+  =  diag(U,U,U  +  C<,0) 
•\nd  A”  =  diag(0,0,0,U  -  C(}. 

4.1.2  Splitting  Based  on  Pressure  Gradient 

The  second  splitting  under  consideration  is  based  on  the  suggestion  by  Vignercn 
(26],  who  split  the  streamwise  preesure  gradient  into  two  parts.  His  discussion  of 
dp/dx  can  also  be  interpreted  as  a  flux  splitting  procedure  given  by, 


pU 

0 

puU  +  u )^xp 
pvU  +  u£vp 

i 

II 

^l*s 

(1  -  w){zp 
(1  -  u){vp 

.  [e  +  p)U  . 

0 

This  splitting  recognizes  that  the  streamwise  ellipticity  arises  from  the  presst 
dic;it  term  inside  the  subsonic  portion  of  the  boundary  layer.  Due  to  this  pressure 
gradient,  downstream  information  can  propagate  upstream.  Thus,  Vigneron’s  no¬ 
tion  was  to  separate  the  streamwise  pressure  gradient  into  parts,  then,  by  choosing 
the  coefficient  w  properly,  one  can  place  those  parts  responsible  for  upstream  influ¬ 
ence  i.ito  the  flux  vector  E~ .  The  vector  E+  therefore  contains  only  those  parts 
governing  the  propagation  from  upstream  to  downstream. 

The  eigenvalues  of  the  Jacobian  matrix  of  E+  in  Eq.  (4.7)  re  the  four  roots 
of  the  following  polynomial, 


(A  -  U)2  { A2  ~  h  +  1  -  w(q r  -  1 ) j t/ A  +  (7  +  u/  -  un)6'2  -  u/C'|)  =  0.  (4,8) 
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The  eigenvalues  of  AJ  are  the  roots  of 

* 

A3(A  -  (1  -7)(1  -u)U\  =0. 


These  roots  are  found  to  be 

A+  =  diag(U,U,\/2{[(~t  +  1)  -  ^  -  1)\U  ±  1  )7U'+AuC\)) 

A"  =  diag(0,0,Q,  -(^  -  1)(1  -  u)U) 

(4.0) 

In  keeping  with  our  purpose,  all  four  eigenvalues  of  A?  =  dE* /dQ  must  be 
positive,  and  those  eigenvalues  of  A/~  =  dE~ /dQ  must  be  negative.  For  supersonic 
flows  ( U  >  C^),  if  we  set  u  equal  to  unity,  then  A+  becomes  A+  =  diag{U,U,U  + 
C(,U  -  Ci)  and  A"  is  identically  zero.  All  entries  of  A+  are  positive,  thus  the 
splitting  is  completed  by  setting  w  =  1  for  supersonic  flows. 

For  subsonic  flows  ( U  <  C^),  the  eigenvalues  of  A,-  are  negative  as  long  as 
u>  <  1.  On  the  other  hand,  the  eigenvalues  of  A?  are  positive  if  the  three  inequalities 

U  >0 

(if  +  1  -  w(-y  -  l)]f/  >  0  .  (4.10) 

(-y  +  w  -  w'yjf/2  -  u/C2  >  0 

are  satisfied.  The  first  inequality  is  straightforward  while  the  second  is  "quivalent 


to 


u  < 


*y  +  i 
^  -  1 


which  is  valid  if  ui  <  1.  The  third  inequality  in  Eq.  (4.10)  results  in 

1 

u<  I  +  (-,  -  l)M» 


(111) 


where  M(  is  the  streamwiso  Mach  number  defined  by  A/<  -■  U/Ct. 

In  summary,  pressure  gradient  splitting  based  on  Vigeneron’s  approach  gives 
the  proper  signs  for  eigenvalues  of  both  Af  and  when  the  following  conditions 


are  met: 
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1.  The  streamwise  contravariant  velocity  it  positive. 

2.  u  it  unity  if  the  flow  is  supersonic. 

3.  u>  satisfies  Eq.  (4.11)  if  the  flow  is  subsonic. 

Three  observations  are  noted  here.  First,  the  flux  vector  splitting  in  Eq.  (4.7)  is 
analogous  to  traditional  Vigneron*based  PNS  procedures.  Second,  the  derivations 
of  the  conditions  for  proper  pressure  gradient  splitting  given  above  differ  from  the 
derivations  given  by  previous  investigators  (3,26),  in  which  the  same  splitting  criteria 
are  obtained  by  considering  the  steady  state  TLNS  equations  in  both  the  inviscid 
and  viscous  limits  of  the  corresponding  simplified  equations.  By  requiring  these 
simplified  equations  to  be  hyperbolic  along  the  streamwise  direction  in  the  inviscid 
limit,  and  to  be  parabolic  in  the  viscous  limit,  they  arrived  at  the  same  conclusions 
given  above.  In  the  present  approach,  the  unsteady  TLNS  equations  are  considered, 
the  splitting  criteria  are  then  obtained  by  taking  into  account  only  the  signs  of 
eigenvalues  of  the  Jacobian  matrix  in  the  streamwise  direction.  Since  the  unsteady 
version  of  the  TLNS  equations  is  hyperbolic  in  time,  by  forcing  these  eigenvalues 
to  be  positive,  we  can  easily  complete  the  splitting  without  considering  the  TLNS 
equations  in  the  inviscid  or  the  viscous  limits  separately.  Third,  in  the  formulation 
of  the  pressure  gradient  splitting,  the  flux  vectors  E*  and  E~  are  defined  without 
a  prior  specification  of  the  matrices  A+  and  A~.  In  fact,  the  flux  vectors  E  +  and 
E~  defined  by  Eq.  (4.7)  are  not  homogeneous;  thus  no  explicit  representations  of 
A +  and  A~  for  the  pressure  gradient  splitting  exist  (as  that  given  in  Eq.  (4.6)). 
However,  if  we  neglect  the  variation  of  u>  with  respect  to  Q,  the  relations 

E+  =  A+Q,  E~  =  ATQ 

similar  to  Eq.  (4.6)  can  be  obtained,  in  which  A*  -  DE*  ;  OQ  and  /t,~  =  0E~  jdQ. 

Just  as  in  the  discussion  of  characteristics  splitting  procedure  in  Section  4.1.1, 
we  have  split  the  flux  vector  into  positive  and  negative  parts. 
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4.2  Obtaining  the  PNS  Procedure  from  the  Navter».Sio.k»s.  Algorithm 


In  the  last  section,  we  have  formulated  general  dux  vector  splitting  for  the 
TLNS  equations.  In  the  special  case  if  the  splitting  is  based  on  characteristics, 
we  have  shown  in  Chapter  3  that  the  resulting  TLNS  equations  can  be  efficiently 
solved  by  approximate  factorization  procedures.  The  traditional  Vigeneron’s  parab¬ 
olization  procedure  has  also  been  interpreted  as  a  special  case  of  the  generalized 
flux-vector  splitting  TLNS  equations.  Based  on  this  interpretation,  the  Vigen¬ 
eron’s  Parabolized  procedure  is  equivalent  to  ignoring  the  reverse  sweep  in  a  specific 
flux-vector  splitting  TLNS  procedure.  This  suggests  that  a  general  parabolization 
method  can  be  devised  based  on  an  arbitrary  flux  splitting. 

This  generalized  parabolization  procedure  can  be  achieved  by  simply  neglecting 
the  parts  of  the  flux  vector  E  contributing  to  upstream  propagation.  If  the  vector 
E~  is  identically  zero  (as  it  is  in  supersonic  flows)  the  algorithm  given  in  Eq.  (4.4) 
describes  an  alternating  procedure  in  one  direction.  For  those  cases  where  E~  Is  not 
zero,  v.ecan  likewise  obtain  a  “marching’’  procedure  by  ignoring  the  contribution  of 
E~ .  In  other  words,  the  streamwise  ellipticity  is  suppressed  by  ignoring  the  elliptic 
parts  of  governing  equations,  thus  the  new  approximate  equations  become  parabolic 
in  the  streamwise  direction.  Again,  we  note  that  we  must  maintain  consistency  on 
both  sides  of  Eq.  (4.4),  so  we  also  drop  the  operator  on  the  left-hand 

side.  With  this  approximation,  the  left  hand  side  of  Eq.  (4.4)  becomes  a  parabolic 
operator, 


{I  -  AID  +  At[—At  +  —  D  -  —(Ih—  Bvl  +  R2  —  D{2)\}AQ  - AiR "  (4.12) 

at,  or i  Or)  dr)  drj 


where  the  modified  residual  R"  is  also  parabolized, 


(4.13) 
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Equation  (4.1?)  defines  a  general  parabolized  procedure  based  upon  any  flux 
splitting  for  which  the  eigenvalues  of  Af  are  positive.  This  implies  that  there  are 
an  infinite  number  of  ways  to  accomplish  the  Parabolized  Navier-Stokes  procedures 
and  Vigneron’o  pressure  gradient  method  is  only  a  special  case  of  these  parabolized 
procedures.  With  this  general  form  of  parabolized  procedures,  the  splitting  based 
upon  characteristics  seems  to  be  more  “natural”  than  the  pressure  gradient  splitting 
in  the  physical  sense.  These  two  special  cases  of  general  parabolized  procedures  are 
considered  in  the  next  sections. 


4.2.1  Pressure  Gradient  Splitting 


In  the  special  case  when  E+  is  given  by  Eq.  (4.7),  Eq.  (4.12)  becomes  the 
traditional  PNS  operator  as  given  by  numerous  authors  (for  example,  Refs. 26*29) 
except  that  the  time  derivative  is  included.  These  time-iterative  PNS  equations 
are  to  be  solved  by  iterations  at  each  streamwise  location.  In  other  words,  since 
Eq.  (4.12)  is  now  a  marching  equation  (this  implies  no  upstream  influences  exist), 
it  is  clearly  better  to  iterate  to  convergence  in  time  at  each  line  before  advancing  to 
the  next  streamwise  station.  To  define  this  time-iterative  procedure  more  precisely, 
Eq.  (4.12)  is  rearranged  as, 


U-MD  +  (1  +  +  ail ±B  -  £(*.£*.,  +  JIJAO  = 


+! %  -  *  -  ' 


(4.1-1) 


where  k  is  0  for  first-order  upwind  differencing  in  the  streamwise  derivative  and  is  1 
for  second  order  differencing.  The  superscripts  *  denote  that  these  quantities  are  to 
he  evaluated  based  on  the  converged  solutions.  The  value  of  AQ  on  the  left-hand 
side  of  Eq.  (4.14)  is  driven  to  machine  accuracy  by  time  marching  at  one  station, 
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and  then  the  procedure  marches  to  the  next  ^-station  and  so  forth.  As  will  be 
demonstrated  later,  this  iteration  can  be  driven  to  machine  accuracy  in  less  than 
10  iterations  for  two-dimensional  problems. 


4.2.2  Characteristics  Splitting 

If  the  flux  vector  J?+  in  Eq.  (4.14)  is  chosen  as  that  given  in  Eq.  (4.6),  a 
similar  time-dependent  PNS  procedure  can  also  be  obtained.  This  marching  pro¬ 
cedure  differs  from  the  more  traditional  pressure  gradient-split  PNS  procedure  in 
the  parabolized  approximation.  As  we  have  seen,  the  pressure  gradient  splitting 
algorithm  omits  parts  of  the  pressure  gradient,  while  the  characteristic  splitting  al¬ 
gorithm  neglects  those  parts  with  upstream-p'opagating  acoustic  wave.  The  latter 
is  more  appropriately  described  by  an  appeal  to  the  physics  of  the  flow.  The  differ¬ 
ences  between  these  two  algorithms  are  also  indicated  by  the  different  eigenvalues 
of  the  Jacobians  of  E  + .  As  will  be  shown  later,  the  calculations  based  on  this  PNS 
procedure  give  results  that  are  almost  identical  to  or  even  slightly  better  than  those 
based  on  the  pressure  gradient  splitting  that  is  traditionally  used. 

4-2.3  Non-Iterative  PNS  Procedure 


The  PNS  algorithms  discussed  above  include  the  temporal  derivative,  while 
in  the  traditional  PNS  procedure,  the  solutions  are  obtained  by  a  simple  space 
marching  without  iterations.  To  obtain  this  marching  procedure,  we  first  re-write 
Eq.  (4.12)  without,  using  the  delta  form.  By  cancelling  terms  on  the  lefl-hnnd 
side  with  those  on  the  right-hand  side  (given  in  Eq.  (4.13)),  the  time-dependenl 
algorithm  becomes 


(«+ +  gf  -  ">  - 


n  +  I 


=  Qn.  (4.15) 


Note  that  for  characteristic  splitting,  does  not  cancel  with  Af  Q  since,  A*  f  A* 
if  the  flow  is  subsonic.  Therefore,  Eq.  (4.15)  is  only  approximately  valid  for  the 
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characteristic  splitting  algorithm,  but  this  is  not  severe  since  in  general  the  subsonic 
layer  is  very  thin.  If  we  allow  At  to  go  to  infinity,  and  use  the  chain  rule, 


dE+  _  dE+  dQ  _  dQ 
d(  dQ  df  ~  A'  at ' 


H  16) 


Equation  (4.15)  then  becomes 


0E+ 

di 


.dQ  dF  -  d  d  d 

=  AlTl=-T,+H  +  rtR'ir,B"+R'T,B")Q- 


(■1.17) 


We  can  now  linearize  each  term  in  Eq.  (4.17)  according  to  the  Taylor  series  expan¬ 
sions, 

Ft+i  =  Ft  +  BAQ 

(418) 

//.+  ,  =  //,  +  DAQ 

where  subscripts  represent  the  £  direction  grid  number,  B  and  D  are  Jacobians  of 
F  and  H ,  and  AQ  is  now  interpreted  in  a  spatial  rather  than  a  temporal  sense, 

a  q  =  gl+1  -Q,. 

With  the  substitution  of  Eq.  (4.18)  into  Eq.  (4.17),  we  have, 


(y|*  -  A«/>+A«|  — fl-  +  R,~B„)\).AQ  = 

_  Afl—  -  H  -  —IR  +  R  dAls\ 

£lclrj  H  a^R'  dr,  R 1  dr,  ' 


(1.19) 


This  equation  now  can  be  used  to  solve  Q,+  |  without  iterations  at  i  +  1  station 
provided  that  Qt  is  given.  Equation  (4.19)  is  referred  to  in  the  literature  as  a 
space-marching  PNS  algorithm.  The  formulation  above  shows  that  any  fiux-vcctor 
splitting  defined  by  Eq.  (4.1)  can  be  used  to  obtain  a  distinct  non-iterative  PNS 
algorithm.  In  particular,  the  pressure  gradient  splitting  noted  above  gives  the  tra¬ 
ditional  PNS  procedure,  Equation  (4.17)  also  chows  that  the  characteristic  splitting 
suggested  by  Steger  and  Warming  can  be  used  to  formulate  a  parabolized  algorithm 
as  well. 
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4.2.4  Comparisons  of  Time-Iterative  and  Space-Marching  Algorithms 

In  the  discussion  above,  both  the  time-iterative  and  space-marching  PNS  pro¬ 
cedures  are  formulated.  By  intuition,  one  might  expect  that  the  space-marching 
algorithm  is  more  efficient  than  the  time-iterative  algorithm  because  it  does  not 
require  local  iterations  at  each  £  station.  There  are,  however,  other  issues  involved. 
To  demonstrate  this,  we  compare  the  differences  between  two  approaches.  First,  we 
consider  the  difference  in  the  final  converged  solutions  of  the  two  methods.  From 
Eq.  (4.12),  if  AQ  is  driven  to  zero,  the  steady  state  solutions  of  the  time-iterative 
procedure  can  be  obtained  as, 


dF 

i-AA+Q)  +  ir-vT-  =  o 


dt 


dr) 


(4.20) 


where  A*Q  =  E+  and  V.T.  is  used  to  represent  the  viscous  terms.  On  the  other 
hand,  the  solutions  of  the  space-marching  algorithm  are, 


A+ 


dQ 

d< 


+ 


(4-21) 


The  r;  derivatives  and  the  viscous  terms  are  exactly  the  same  for  both  methods.  The 
difference  lies  in  the  (  derivative,  in  which  the  time-iterative  method  utilizes  the 
conservative  form,  while  the  space-marching  method  employs  a  non-conservative 
form. 

To  explore  this  difference  further,  we  compare  the  finite-difference  representa¬ 
tions  of  the  {  derivative  in  Eqs.  (4.20)  and  (4.21).  For  clarity,  we  restrict  to  first 
order  in  £.  Upon  discretization,  the  £  derivative  in  Eq.  (4.20)  becomes 

+  i  ~  M  +  Q), 


«-<?.- 


Qt 


while  the  one  in  Eq.  (4.21)  is 
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Because  the  first  row  of  the  matrix  A+  is  related  to  the  continuity  equation,  the  non- 
conservative  form  in  Eq.  (4.21)  can  be  expected  to  give  mass  conservation  errors. 
For  a  grid  which  has  no  stretching  in  the  £  direction,  these  mass  errors  are  not  severe 
since  the  metric  coefficients  are  constants  in  the  £  direction.  However,  when  grid 
stretching  is  used  in  £,  the  mass  error  can  be  expected  to  accumulate  with  (  because 
of  the  variation  in  the  metrics.  Numerical  experiments  using  the  non-iterative 
(space-marching)  scheme  show  that  a  global  mass  error  of  order  one  is  observed  for 
even  a  moderately  stretched  gr»  J.  For  highly  stretched  grids,  this  accumulation  leads 
to  numerical  instability.  On  the  other  hand,  the  time-iterative  algorithm  worked 
well  with  either  uniform  or  highly  non-uniform  grids.  Consequently,  the  penalty 
paid  for  local  iterations  in  the  time-iterative  algorithm  can  be  at  least  partly  offset 
by  using  a  stretched  grid. 

An  alternative  procedure  for  ensuring  mass  conservation  with  non-iterative 
scheme  on  stretched  grids  has  been  proposed  by  Schiff  and  Steger  |43|,  although  it 
does  not  appear  to  have  been  widely  used.  In  their  approach,  instead  of  directly 
using  the  chain  rule  given  by  Eq.  (4.16)  to  represent  dE+  /d£  in  Eq.  (4.17),  the  flux 
vector  E  is  linearized  before  discretizing.  This  procedure  is  obtained  by  noting  that 
the  flux  vectors  at  two  consecutive  locations  can  be  linearized  according  to 


EZt=E?+AtiQt+l-Qt) 


By  using  Eq.  (4.22)  and  Eq.  (4.23),  we  have 

d£  "  AC'4,  *Z{A' 


(1.22) 

(1.23) 


(1.24) 


With  the  use  of  Eq.  (2.24)  in  Eq.  (4.17),  the  space-marching  PNS  algorithm  gives 
better  mass  conservation.  Numerical  experiments  with  this  approach  prove  to  he 
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able  to  conserve  mass  within  an  error  of  1%  for  a  non-uniform  grid  in  a  moderate 
expansion  ratio  (around  30)  nozzle.  However,  for  more  realistic  problems  such  as 
flows  through  the  272  :  1  nozzle  investigated  in  Chapter  3,  the  75  axial  grid  lines  had 
to  be  increased  to  300  to  enable  the  modified  space-marching  algorithm  to  match 
the  conservative  time-iterative  method  with  75  axial  grid  lines  in  accuracy.  The 
space-marching  procedure  without  the  SchifT-Steger  modification  led  to  global  mass 
errors  of  more  than  50%  even  with  300  axial  grid  lines. 

The  second  difference  to  be  addressed  is  the  requirement  of  a  safety  factor, 
d,  in  defining  the  parabolized  operator.  Parabolized  Navier-Stokes  calculations 
reported  in  the  literature  (for  example,  Ref.  (29|)  traditionally  use  a  safety  factor 
in  Eq.  (4.11).  This  results  in 

U<  l  + 

where  the  safety  factor  9  is  generally  chosen  as  0.85  or  smaller.  Numerical  experi¬ 
ments  with  the  space-marching  algorithm  indicate  that  9  can  not  be  greater  than 
0.85  without  numerical  instability.  On  the  other  hand,  with  the  use  of  the  time- 
iterative  algorithm,  9  can  always  be  set  equal  to  unity.  The  results  presented  in  the 
next  section  also  show  that  the  solutions  with  9  =  1  are  more  accurate  than  those 
with  9  -=  0.85  as  compared  to  the  Navier-Stokes  solutions.  From  the  derivation  of 
u>  given  in  Section  4.1.2,  it  is  clear  that  there  is  no  theoretical  reason  for  requiring  a 
safety  factor.  The  necessity  of  a  safety  factor  in  the  space-marching  method  is  only 
to  make  the  algorithm  stable,  and  is  not  inherent  with  the  parabolized  equations. 

As  a  final  comparison,  we  note  that  the  left  hand  side  operator  of  the  time- 
iterative  procedure  is  more  diagonally  dominant  than  that  of  the  space-marching 
algorithm.  In  fact,  vanishing  elements  appear  on  the  diagonal  of  the  left  hand 
side  matrix  in  the  space-marching  algorithm  due  to  the  absence  of  the  identity 
matrix  /  in  Eq.  (4.19).  As  a  consequence,  pivoting  strategies  are  required  to  solve 
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Eq.  (4.10).  Contrarily,  the  diagonal  elementa  are  always  non-zero  for  the  lime- 
iterative  algorithm,  and  pivoting  is  not  required.  Our  experience  shows  that  a  30% 
saving  in  computational  time  per  iteration  is  gained  by  the  solution  of  a  block  tri¬ 
diagonal  matrix  without  pivoting.  Again,  this  difference  would  help  to  make  the 
time-iterative  PNS  procedure  more  economically  competitive  with  the  non-iterative 
procedure. 

4.2.6  Stability  Analysis  of  Time-Iterative  PNS  Algorithms 

To  validate  the  time-iterative  algorithm  developed  above,  the  Fourier  stability 
analysis  of  Eq.  (4.14)  is  given  as  follows.  The  amplification  matrix  of  the  variable 
Q  is  defined  •' 

Qn+l  =  GQn. 

From  a  von  Neumann  analysis,  G  can  be  found  to  be 

L\G  =  / 

where  the  matrix  L\  is, 

Af  .  Af  Af 

L\  =  /  -  Af£>  +  ——A,  +  i  —  B  a\nun  +  2- — ~[R\BV\  +  /?2B„2)(1  -  coswJ 
A£  A  n  A  r?2 

and  u ;n  is  the  rj  direction  wavenumber.  Figure  46  shows  maximum  eigenvalues  of 
the  G  matrix  versus  wavenumber  for  typical  supersonic  and  subsonic  conditions. 
The  results  show  that  Eq.  (4.14)  is  unconditionally  stable,  and  that  rapid  conver¬ 
gence  can  be  expected  for  high  values  of  CFL.  The  stability  results  given  above 
are  for  the  pressure  gradient  splitting  method.  Those  of  characteristic  splitting  are 
qualitatively  the  same. 

4.2.6  Results  and  Discussion 


We  have  discussed  two  parabolized  procedures  so  far.  Now,  the  question  to  ask 
is  which  method  is  better  in  terms  of  both  computational  efficiency  and  solution 


Maximum  eigenvalue  of  G 
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accuracy?  To  answer  this  question,  a  series  of  numerical  computations  were  done 
and  their  solutions  were  compared  to  those  from  the  TLNS  algorithms  provided 
in  Chapter  3.  Before  showing  the  results  of  these  comparisons,  we  note  here  that 
all  PNS  algorithms  given  above  require  boundary  conditions  in  the  cross-stream 
direction.  The  procedures  at  the  wall  and  the  centerline  discussed  in  Section  3.2.1 
are  equally  applicable  to  PNS  algorithms.  At  the  starting  plane,  the  flow  variable 
Q  must  be  given. 

The  test  problem  for  the  comparisons  following  is  again  the  (low  through  the 
high  expansion  ratio  nozzle  with  an  area  ratio  of  272  :  1.  The  same  flow  conditions 
and  properties  described  in  Section  3.6  were  used,  including  the  inlet  conditions,  the 
ratio  or  specific  heats,  and  the  75  x  50  grid  system  (for  the  non-iterative  algorithm, 
a  more  refined  300  x  50  grid  was  used).  The  Reynolds  number  was  taken  to  be 
1.4  x  10*  based  on  the  throat  radius  for  all  calculations  presented.  All  flowfield 
results  presented  are  for  laminar  calculations. 

The  numerical  efficiency  of  the  time-iterative  PNS  procedure  is  shown  in  Fig.  47 
for  representative  conditions.  This  figure  shows  the  convergence  at  a  specific  £  sta¬ 
tion  by  plotting  the  L-2  norm  of  uQ/Q  associated  with  the  four  equations  (con¬ 
tinuity,  momentum,  and  energy  equations)  as  a  function  of  the  iteration  number. 
Both  inviscid  and  viscous  results  are  shown  on  Fig.  47.  A  CFL  number  of  10°  was 
chosen  for  both  cases.  The  convergence  clearly  indicates  that  machine  accuracy 
was  reached  in  less  than  10  iterations  and  the  inviscid  case  converges  slightly  faster 
than  the  viscous  case.  Acceptable  convergence  (corresponding  to  a  reduction  of 
six  order  of  magnitude  in  the  L-2  norm)  is  obtained  in  4  iterations.  As  we  can  see 
from  Eq.  (4.14),  when  the  time  step  goes  to  infinity,  the  time-iterative  algorithm 
approaches  Newton’s  method.  Hence,  the  quadratic  convergence  shown  in  Fig.  47 
is  to  be  expected.  This  rapid  convergence  has  been  generally  observed  for  all  PNS 


L-2  Norm  of  A Q/Q 


Number  of  iterations 


Figure  47.  Convergence  rate  of  time-iterative  PNS  procedure  at  typic'  I 
axial  location 
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calculations  to  date.  The  convergence  shown  in  Fig.  47  is  based  upon  pressure 
gradient  splitting,  but  is  also  representative  for  characteristic  splitting  algorithm. 

The  solutions  of  PNS  algorithms  are  compared  to  those  from  thin-layer  Navier- 
Stokes  calculations  in  Fig.  48.  The  upper  plot  s/'v*  the  Mach  number  contours 
by  using  the  pressure  gradient  splitting  PNS  method,  the  lower  plot  shows  similar 
results  for  the  TLNS  solutions.  This  comparison  indicates  that  the  PNS  procedure 
gives  solutions  that  are  almost  identical  to  those  of  the  TLNS  equations.  As  will 
be  shown  later,  the  characteristic  splitting  also  gives  results  that  are  even  closer  to 
the  TLNS  solutions 

To  further  compare  the  flowfield  details,  the  pressure  distribution  and  stream- 
wise  velocity  profiles  at  the  exit  plane  are  plotted  in  Fig.  49  and  Fig.  50.  Each  figure 
compares  four  different  procedures.  They  are  the  TLNS  method,  the  pressure  gra¬ 
dient  splitting  PNS  with  0  =  1.0  and  d  =  0.85,  and  the  characteristic  splitting 
PNS.  The  pressure  profiles  in  Fig.  49  show  that  the  computed  pressure  distribution 
by  pressure  gradient  splitting  without  safety  factor  is  almost  identical  to  that  by 
characteristic  splitting,  except  the  former  slightly  overshoots  the  pressure  at  the 
centerline.  Both  methods  agree  very  well  with  the  TLNS  results  and  the  charac¬ 
teristics  splitting  method  is  slightly  better  than  the  pressure  gradient  PNS.  Figure 
49  also  shows  that  the  use  of  a  safety  factor  of  0,85  in  pressure  gradient  splitting 
deteriorates  the  solution  accuracy.  As  can  be  seen,  the  use  of  the  safety  factor 
causes  about  a  25%  undershoot  in  the  pressure  at  the  centerline  and  about  a  5% 
overshoot  at  the  wall. 

Similar  comparisons  associated  with  the  velocity  profiles  shown  in  Fig.  50  show 
that  all  four  procedures  predict  fairly  close  velocity  distributions.  Correct  values  of 
velocity  together  with  incorrect  values  of  static  pressure  indicate  that  the  entropy 
(stagnation  pressure)  is  not  well  conserved,  a  phenomenon  frequently  encountered 
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Figure  50.  Comparison  of  cross-stream  velocity  profile  at  exit  plane  for 
various  PNS  results  with  TLNS  calculations 
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in  numerical  schemes. 

The  effect  of  setting  9  =  0.85  is  also  indicated  in  Fig.  51,  which  n 

number  contours  for  this  case.  Comparison  of  this  figure  with  the  t.n.  .usc.ions 
given  in  Fig.  48  shows  that  the  over-suppression  of  the  streamwise  pressure  gradient 
by  using  safety  factor  other  than  one  alters  the  flowfteld  structure  dramatic  ally. 

Further  comparisons  are  shown  in  Fig.  52  and  Fig.  53,  where  the  wall  pressure 
distributions  and  the  locations  of  the  sonic  line  along  the  axial  direction  are  plotted. 
Figure  52  shows  that  all  four  methods  give  almost  identical  wall  pressure  distribu¬ 
tions  but  again  the  case  with  the  safety  factor  included  is  not  quite  as  accurate. 
The  sonic  line  locations  shown  in  Fig.  53  are  obtained  by  interpolation  between  grid 
points.  This  figure  shows  that  the  three  PNS  algorithms  give  basically  the  same 
subsonic  layer  thickness  (the  distance  from  the  wall  to  the  sonic  point).  Although, 
the  O  =  0.85  case  underpredicts  the  thickness  of  subsonic  layer  by  about  1%,  which 
is  the  worst  among  the  three  algorithms.  These  results  indicate  that  the  PNS  ap¬ 
proximation  gives  solutions  that  are  acceptable  in  accuracy  for  the  high  Reynolds 
number  flow  without  separation,  as  in  current  test  problem. 

The  discussion  above  demonstrates  that  for  better  solution  accuracy,  t,»e  safety 
factor  should  not  be  less  than  one  (which  is  easily  done  by  using  the  time-iterative 
algorithm).  Furthermore,  the  characteristic  splitting  PNS  procedure  gives  solutions 
that  are  as  accurate  as,  or  even  more  accurate  than  (as  in  current  test  problem)  the 
traditional  pressure  gradient  splitting  PNS  method. 

So  far,  all  results  shown  for  supersonic  viscous  calculations  including  hot  It  from 
TLNS  and  PNS  algorithms  are  obtained  by  using  second  order  differencing  in  the  £ 
direction.  To  demonstrate  the  difference  in  accuracy  between  first  order  and  second 
order  accurate  upwind  differencing  in  two  dimensions,  the  first-order  PNS  results 
of  the  same  test  problem  (flows  through  the  high  expansion  nozzle)  are  shown  in 


F.gure  51.  Mach  number  contours  computed  by  PNS  algorithm  based  on  pressure 
gradient  splitting  and  using  a  safety  factor  of  0.85 


ressure  distribution  on  the  wall  for  TLNS  and  PNS  procedures 


subsonic  layer  as  a  function  of  axial  distance  along  the 
vine  results  of  TLNS  and  various  PNS  calculations 


142 


Fig.  54.  Comparisons  of  Fig.  48  and  Fig.  54  show  that  the  oblique  shock  wave  from 
the  first-order  solutions  is  not  as  sharp  as  that  from  second-order  solutions  due  to 
the  smearing  effect  resulting  from  the  inherent  second  order  artificial  dissipation  of 
the  first-order  upwind  differencing.  Therefore,  for  better  solution  accuracy,  second 
order  upwind  differencing  should  always  be  used. 

4.3  Global  PNS  Procedures 

For  flows  with  strong  upstream  influences  such  as  separated  flows,  the  marching 
type  PNS  procedure  as  discussed  in  Section  4.2  can  lead  to  serious  errors  in  the 
numerical  solution  due  to  the  suppression  of  the  streamwise  ellipticity.  To  allow  the 
upstream  propagation  of  acoustic  wave  inside  the  subsonic  layer,  thus  preserving 
the  streamwise  elliptic  behavior,  the  dE~ /d£  term  in  Eq.  (4.2)  must  be  included; 
thus  the  Navier-Stokes  procedures  discussed  in  Chapter  3  must  be  used  instead  of 
the  parabolized  algorithms  provided  in  this  chapter. 

In  the  traditional  PNS  approach,  numerous  attempts  have  been  made  to  take 
into  account  the  upstream  influences  by  identifying  global  pressure  iterations.  The 
basic  idea  of  global  pressure  iterations  is  to  update  the  pressure  field  by  providing 
some  sort  of  stable  differencing  scheme  for  the  the  omitted  (1  -u>)dp/d£  term.  This 
is  usually  done  by  evaluating  (1  —u)dp/d(  from  a  forward  difference  and  using  the 
updated  value  of  pressures  at  downstream  locations,  as  suggested  in  the  works  of 
Rakiclt  |44),  and  Lin  and  Rubin  (45).  Davis  and  co-workers  (46)  and  Darnel  t  and 
Davis  (47)  also  developed  a  global  pressure  iteration  procedure  by  appending  a 
fictitious  unsteady  term,  dp/ dr,  on  the  steady  state  equation,  then  updating  the 
pressure  field  by  a  two-step  alternating  direction  explicit  (ADE)  procedure.  These 
global  pressure  iterations  are  summarized  by  Thompson  and  Anderson  (48|. 

In  the  present  study,  by  obtaining  the  PNS  procedure  from  the  TLNS  equa¬ 
tions,  the  procedure  for  incorporating  a  global  pressure  iteration  procedure  becomes 


r 
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Figure  54.  Mach  number  contours  computed  by  PNS  algorithm  based  on  pressure 
gradient  splitting,  first-order  upwind  results. 
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obvious.  We  need  only  return  to  the  complete  TLNS  equations  (Eq.  (4.2)).  There¬ 
fore,  all  thr'ee  approximate  factorization  algorithms  provided  in  Chapter  3  can  be 
interpreted  as  global  pressure  iteration  procedures. 

The  TLNS  algorithms  (or  so-called  global  pressure  iteration  procedures  in  the 
traditional  PNS  approach)  developed  in  this  study  are  based  upon  the  approximate 
factorization  of  the  TLNS  equations,  and  therefore  have  both  physical  and  mathe¬ 
matical  connections  to  the  equations  of  motions,  while  the  global  pressure  iteration 
algorithms  are  concerned  with  arbitrary  iterative  processes  for  the  pressure  gradi¬ 
ent  (1  -  u)dp/d£,  which  are  unrelated  to  the  physical  equations,  as  a  consequence, 
some  sort  of  relaxation  scheme  is  required.  This  suggests  that  the  mathematically 
and  physically  well-behaved  TLNS  algorithms  based  on  approximate  factorization 
can  be  used  instead  of  global  pressure  iteration  procedures. 

As  an  example  of  the  interpretation  of  global  pressure  iterations  based  on  TLNS 
algorithms,  the  following  procedure  is  suggested: 

1.  Obtain  an  initial  PNS  solution  by  marching  from  upstream  to  downstream 
using  Eq.  (4.4). 

2.  Solve  the  discretized  equation  of  the  PNS-AD!  algorithm,  Eq.  (3.40),  by  the 
following  two  equations, 

{/  -  At D  +  At\~A?  +  j^B-  -^(Ri  j^BvX  +  /?2~Bv2)|}AQ’  =  -A tfi' 

(4.25) 

(/  -  A  ID  +  Af^ADAQ  =  (/  -  AfZ?)AQ'  (4.26) 

3.  Update  the  dependent  variable  Q,  according  to, 

Qn+ 1  =  Qn  +  AQ 


until  the  converged  steady  state  is  reached. 
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The  first  step  is  used  to  obtain  a  good  initial  condition  for  the  TLNS  calculations. 
The  first  equation  in  the  second  step,  Eq.  (4.25),  is  equivalent  to  the  time-iterative 
PNS  algorithm  and  the  second  equation,  Eq.  (4.26),  is  augmented  in  order  to  provide 
a  mechanism  to  allow  upstream  propagation  to  take  place  inside  the  subsonic  layer. 
In  the  supersonic  region,  A is  identically  2ero  and  the  left  hand  side  operator  in 
Eq.  (4.26)  reduces  to  an  identity  matrix;  hence,  only  Eq.  (4.25)  has  to  be  solved. 

Figure  55  shows  typical  convergence  of  the  TLNS  procedure  mentioned  above 
when  applied  to  the  high-expansion  ratio  nozzle  calculation  given  in  Section  4.2.6. 
It  requires  only  110  iterations  to  reach  machine  accuracy;  acceptable  convergence 
is  achieved  in  25  iterations. 


CHAPTER  E 


THE  APPLICATION  OF  TIME-ITERATIVE  SCHEMES  TO 
VISCOUS  SWIRLING  NOZZLE  FLOWS 


As  examples  of  the  application  of  the  Nr  ier-Stokes  algorithms,  swir'ing  vis¬ 
cous  flows  in  transonic  and  supersonic  propulsion  nozzles  are  investigated  in  this 
char'er.  The  central-differenced  ADI  and  the  flux- vector  splitting  algorithms  dis¬ 
cussed  in  Chapter  3  are  utilized  to  solve  the  thin-layer  Navier-Stokes  equations  for 
axisymmetric  two-dimensional  flows  wit'  swirl.  The  effects  of  swirl  on  viscous  flows 
are  identified  for  nozzles  with  mild  to  high  expansion  ratios.  Both  flowfield  details 
and  the  integral  nozzle  performance  are  compared  to  previously  publ  siied  inviscid 
calculations. 

G.l  Governing  Equations  and  Numerical  Algorithms 


The  swirling  nozzle  flow  inside  an  axisymmetric  nozzle  can  be  described  by  the 
three-dimensional  Navier-Stokes  equations.  If  we  assume  the  flow  is  axisymmetric. 
all  circumferential  derivatives  vanish  and  the  system  of  equations  reduces  to  two 
dimensions.  The  resulting  Navier-Stokes  equations  in  vector  form  can  be  written  as 


dQ  dE  dF  rr  dEv  dFv 

—  +  —  +  —--H  +  ~r—  +  -r— . 
dt  dx  dy  dx  dy 


(5-1) 


where  x  and  y  are  the  axial  and  radial  coordinates,  respectively.  The  flow  variab'e 
Q  in  Eq.  (5.1)  is  defined  as 
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in  which  u,  t>,  and  w  represent  axial,  radial,  and  circumferential  velocity  compo¬ 
nents.  Th*  'inviscid  flux  vectors  E  and  F  are  given  by 


+  P 

E  —  y  puv 
puw 

.(e  +p)u. 


pv 

puv 

F  =  y  pv7  +  p 
pvw 

-(*  +p)u. 


Vpcoua  terms  are  included  in  flux  vectors  Ev  and  Fv% 


Ev  =  y  + 

.Mv(fe  +  fc)  +Mu(|ft  -  !fe)  +  *Jf . 


Fv  =  y 


Mils- 3S5) 


U-cfc  +  is) + -  its) + + *k 


The  source  vector  //  is  defined  by 


-!*M 

H=  p  +  pw2  -  fpj  +  -  3t'f? 

—pvw  ~ 

.-3^(^uv)  -  I  ^(mv2)  -  2Aiw|5  -  u>2§*. 

The  system  of  equations,  Eq.  (5.1),  is  similar  to  the  Navier-Stokes  equations  in 
axisymmetric  two  dimensions  (Eq.  (3.2))  except  the  tangential  momentum  equal  ion 
is  included  in  Eq.  (5.1)  to  take  into  account  the  variation  of  the  cirr umferenlinl 
velocity  in  axial  and  radial  directions. 

Following  similar  procedures  given  in  Section  3.1,  the  thin-layer  version  of 
Eq.  (5.1)  in  general  coordinate  system  for  two-dimensional  swirling  flows  can  be 


written  as 


dQ  +  dE_  dF_  _  -  dT\ 

dt  ^  drj  +  dr) 


(5.2) 
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where  the  flow  variable  and  the  flux  vectors  are 


-  y  j, 

Q  =  j(p,pu,pv,pw,c)  , 


and 


'  pU  * 

■  pV 

7 

puu  +  itp 
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-  du  ,  du 
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Q  3  «  ,  Sjil  4.  2ii£l  §*1  +  i.^1  4.  Q,  3uu 

L  Q<  Tj  p  ^  2  TfT  ^  2  IT  ^  2  drj  ^  2  a»»  J 

In  which  <*i  through  a4  follow  the  same  definitions  described  in  Chapter  3  and  ag 
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<*6  =  p{ril  +  f?J). 

The  source  vector  H  for  swirling  flows  is  simply  H  =  H/J.  Again,  if  p  is  set  equal 
to  zero  and  dFv/dq  is  omitted,  Eq.  (5.2)  reduces  to  the  Euler  equations  which 
describe  inviscid  swirling  flows. 

Equation  (5.2)  takes  the  same  form  as  Eq.  (3.9)  except  for  the  additional  en¬ 
tries  arising  from  the  tangential  momentum  equation.  Therefore,  all  numerical 
algorithms  discussed  in  Chapter  3  and  Chapter  4  are  presumably  applicable  for 
the  present  governing  equations.  According  to  the  nature  of  the  flow,  different  nu¬ 
merical  algorithms  will  be  employed  to  solve  transonic  and  supersonic  flows.  For 
transonic  flows,  we  choose  the  implicit  ADI  procedure  instead  of  MacCormack’s 
explicit  algorithm  which  was  used  by  previous  workers  (32-34 ] .  The  details  of  this 
implicit  ADI  ocherne  have  been  discussed  in  Section  3.2  and  will  not  be  repeated 
here. 
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As  mentioned  earlier,  the  implicit  ADI  procedure  becomes  inefficient  and  some¬ 
times  even  leads  to  numerical  instability  if  a  large  portion  of  the  flowfield  is  super¬ 
sonic  For  this  viscous  supersonic  swirling  flow,  the  hybrid  upwind-central  differ¬ 
encing  algorithms  described  in  Section  3.4  can  be  chosen.  For  swirling  flows  where 
upstream  influence  is  not  significant,  the  parabolized  procedures  discussed  in  Chap¬ 
ter  4  are  also  applicable.  The  major  difference  In  numerical  procedures  between  the 
present  swirling  flow  solvers  and  the  non-swirling  axisymmetric  solvers  discussed  in 
Chapter  3  is  that  the  block  size  of  the  left  hand  side  matrix  for  the  present  case  is 
5x5  while  that  for  the  non-swirling  case  is  4  x  4.  Therefore,  numerical  procedures 
for  the  swirling  flows  are  more  time-consuming  than  those  for  the  non-swirling  cases. 

5.2  Boundary  Conditions 

Previously  defined  boundary  procedures  can  be  easily  extended  to  swirling  flow 
calculations.  For  supersonic  flows  in  the  meridian  plane  at  the  inlet,  the  flow  variable 
Q  i3  completely  specified.  For  subsonic  inflows  at  the  inlet,  the  stagnation  tempera¬ 
ture,  the  stagnation  pressure,  the  meridian  plane  streamline  angle  ifr  =  tan~ 1  ( u / u ) , 
and  the  swirl  angle  <t>  =  tan~l  (w/u)  are  specified,  the  remaining  one  unknown 
comes  from  the  characteristic  equation  corresponding  to  the  single  negative  eigen¬ 
value.  The  swirl  angle  profile  at  the  inlet  is  assumed  to  be  one  of  constant  angle, 
free  vortex,  or  forced  vortex,  which  are  the  same  as  in  Ref.  (34)  except  that  the  swirl 
angle  asymptotically  approaches  zero  at  the  wall  for  all  viscous  calculations. 

At  the  wall,  four  characteristic  equations  and  the  tangency  condition  are  em- 
posed  for  inviscid  calculations,  while  no-slip  conditions  together  with  zero  normal 
pressure  gradient  and  isothermal  or  adiabatic  wall  conditions  are  used.  Symmet¬ 
rical  conditions  are  applied  at  the  centerline.  At  the  exit,  either  extrapolation  or 
fixed  back  pressure  conditions  can  be  chosen  as  described  in  Chapter  3. 


To  give  assessments  of  the  nozzle  performance,  several  integral  parameters  are 
defined  as  the  following.  These  include  the  discharge  coefficient  Co\  the  vacuum 
stream  thrust  efficiency  rjvt ,  the  specific  impulse  efficiency  rjsj ,  and  the  nozzle 
flowfield  as  a  function  of  the  inlet  swirl  number  S,, 

CD  s  -  2  f  puydy/ivlt  -  vlt){p‘u'),d 

^1  d  J  y€ 
f  V*** 

T)vi  =  2  /  (p  +  pu7)ydy/ {y7wt  -  y]t){pt  +  ptu\)xi 


1st  = 


r v rVv 
>.  =  /  puwy7dy/ ywl  /  pu7ydy. 

•'Vt.  Jye, 


The  subscripts  i,  t,  e,  c,  w,  and  id  denote  inlet,  throat,  exit,  centerline  or  centerbody, 
wall,  and  ideal  conditions,  respectively.  The  quantity  m  is  the  mass  flow  rate.  The 
ideal  conditions  are  obtained  from  one-dimensional  isentropic  values  at  the  same 
stagnation  conditions  as  the  actual  flow.  The  discharge  coefficient  can  be  interpreted 
as  a  measure  of  the  loss  in  mass  flow  rate  due  to  two-dimensionality  and  the  swirl 
The  swirl  number  is  defined  as  the  ratio  of  the  axial  flux  of  flow  angular  momentum 
divided  by  the  axial  flux  of  axial  momentum  and  is  a  direct  measure  of  the  level  of 
swirl  at  the  nozzle  inlet. 


5.3  Nozzle  Flowfield 
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To  place  present  viscous  swirling  calculations  in  perspective  with  previously 
published  inviscid  results,  the  implicit  ADI  scheme  is  applied  to  calculate  transonic 
flows  through  a  convergent-divergent  nozzle,  an  annular  plug  nozzle  and  a  converg¬ 
ing  nozzle,  which  all  have  been  investigated  by  Dutton  (34).  As  an  example  of 
predominantly  supersonic  flowfield,  the  viscous  swirling  flow  through  the  272  :  1 
contoured  nozzle  previously  given  is  calculated  by  using  the  upwind-central  differ¬ 
encing  algorithms.  Only  laminar  results  are  shown  for  all  three  transonic  case-,  and 
both  laminar  and  turbulent  results  are  presented  for  the  contoured  nozzle. 


The  35°  -  18.5°  convergent-diverging  nozzle  calculated  by  Dutton  is  analyzed 
in  the  first  series  of  computations.  The  geometry  of  this  nozzle  is  shown  in  Fig.  50 
and  Fig.  57  for  inviscid  and  viscous  computations,  respectively.  The  63  x  30  equally 
spaced  grid  in  both  the  axial  and  radial  directions  as  shown  in  Fig.  56  is  for  inviscid 
calculations,  while  the  63  x  50  grid  with  strong  clustering  near  the  wall  as  shown 
in  Fig.  57  is  for  viscous  calculations. 

The  convergence  rates  of  inviscid  and  viscous  cases  are  shown  in  Fig.  58  and 
Fig.  59,  respectively.  For  inviscid  computations,  as  shown  in  Fig.  58,  the  L-2  norm 
of  AQ/Q  reduces  9  orders  of  magnitude  in  250  iterations  for  the  zero-swirl  case, 
which  is  typical  for  a  ADI  scheme.  Also,  the  presence  of  swirl  is  seen  to  slow  down 
the  convergence  rate  substantially.  The  convergence  rate  for  the  viscous  calculation 
is  dominated  by  the  boundary  layer  near  the  wall,  hence  it.  is  in  general  slower  than 
that  of  inviscid  calculations  as  is  seen  in  Fig.  59  (for  300  iterations,  the  L-2  norm 
drops  only  four  orders  of  magnitude).  These  results  show  that  the  convergence  for 
viscous  calculations  is  insensitive  to  swirl. 

Calculations  of  the  flow  in  the  converging-divergent  nozzle  have  been  completed 
for  a  number  of  nozzle  Reynolds  number  conditions  including  the  inviscid  case.  Fig¬ 
ures  60  and  61  compare  Mach  number  contours  for  the  no-swirl  and  the  swirl  cases 
for  the  inviscid  and  one  of  the  low  Reynolds  number  viscous  calculations,  respec¬ 
tively.  The  inviscid  results  are  in  good  agreement  with  Dutton's  calculations.  The 
viscous  case  is  for  a  Reynolds  number  of  7000  based  on  the  inlet  radius  and  inflow 
properties.  These  viscous  Mach  number  contours  indicate  that  the  introduction  of 
swirl  primarily  affects  the  axial  velocity  near  the  centerline  as  in  the  inviscid  case, 
although  there  are  some  changes  beginning  to  occur  near  the  wall  in  the  diverging 
section. 


Figure  56.  [nvjscid  grid  for  35°-18.5°  convergent-divergent  nozzle 
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The  integral  nozzle  performance  in  the  presence  of  viscosity  is  plotted  in  Fig.  62 
against  the'swirl  number.  This  figure  shows  the  discharge  coefficient  (Cp)  and  the 
vacuum  stream  thrust  efficiency  4  (r?v»)  ns  a  function  of  swirl  numbers  for  free 
vortex,  forced  vortex,  and  constant  angle  inlet  swirl  profiles.  The  predicted  Cp  and 
rjv«  values  are  about  2%  less  than  those  of  inviscid  calculations  of  Dutton  at  the 
Reynolds  number  of  7000.  At  a  given  swirl  number,  the  reduction  in  Cp  and  r?,,,  is 
most  prominent  for  the  free  vortex  case  because  a  relatively  larger  swirl  .igle  must 
be  specified  near  the  centerline  in  order  to  achieve  the  same  swirl  number.  This 
larger  swirl  angle  sresults  in  a  larger  reduction  in  the  mass  flow  rate.  A  similar 
phenomenon  was  also  observed  in  Dutton’s  inviscid  calculations  (34).  The  specific 
impulse  efficiency  [vsi)  for  the  viscous  case  is  essentially  constant  and  is  similar  to 
Dutton’s  inviscid  results,  except  the  value  is  0.965  instead  of  0.971. 

The  effects  of  Reynolds  numbers  are  shown  in  Fig.  63.  The  computed  C p  and 
rjva  values  as  functions  of  the  Reynolds  number  are  plotted  for  5,  =  0  and  5,  = 
0.361.  A  constant  angle  swirl  profile  was  used  for  these  computations.  Asymptotic 
values  obtained  from  present  inviscid  calculations  and  from  Dutton's  calculations 
are  given  on  the  right.  As  the  Reynolds  number  increases,  Cp  and  rjvt  approach 
the  values  of  inviscid  calculations.  These  results  show  the  degree  of  error  incurred 
by  making  the  inviscid  assumption  for  high  Reynolds  number  flows. 

5.3.2  Convergent  Nozzle  and  Plug  Nozzle 

Viscous  calculations  ar:  proceeded  with  the  elliptically  contoured  converging 
nozzle  and  the  annular  plug  nozzle.  A  small  portion  of  the  wall  has  been  appended 
to  the  converging  nozzle  after  the  throat,  such  that  the  flow  at  the  downstream 
boundary  is  predominantly  supersonic  in  the  streamwise  direction.  This  allows 
easier  implementation  of  the  downstream  boundary  conditions  and  does  not  alter 
the  flowfield  before  the  throat.  The  resulting  Mach  number  contours  are  shown  in 


Forced  Vortex 


Dependence  of  integral  performance  on  inlet  swirl  number  for  C-D 


SfO.361 


16J 


Figure  63.  Dependence  of  nozzle  performance  on  Reynolds  number  for  Si 
and  S,  =  .361 
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Fig.  64  for  5,  =  0  and  S,  =  0.43.  Corresponding  performance  curves  are  plotted  in 
Fig.  65.  THe  predicted  Cp  and  r]vt  values  are  about  1%  less  than  those  of  Dutton’s 
inviscid  calculations  over  the  entire  range  of  inlet  swirl  numbers  for  a  Reynolds 
number  of  1.1  x  104.  The  values  of  ns/  are  about  0.7%  less  than  those  of  inviscid 
results. 

The  results  of  the  annular  nozzle  are  shown  in  Fig.  66  for  viscous  calculations 
for  a  Reynolds  number  of  10*.  The  flowfield  of  the  high  swirl  case  5,  =  1.706  is  very 
different  from  that  of  the  zero  swirl  case  near  the  inlet  region  due  to  the  combined 
effect  of  boundary  layer  and  the  inlet  twirling.  The  total  Mach  number  contours 
for  this  highly  swirled  viscous  flow  differ  from  Dutton’s  inviscid  results  due  to  the 
viscous  effect  on  the  circumferential  velocity.  This  discrepancy  demonstrates  the 
importance  cf  v  scous  anoly-is  for  low  Reynolds  number  flows.  About  4%  of  the 
reduction  in  C'c  and  r/vt  compared  to  Vhe  inviscid  case  can  be  observed  in  Fig.  67. 
Again,  the  reduction  in  specific  impulse  efficiency  lor  viscous  calculations  is  less 
than  that  for  the  inviscid  results. 

5.3.3  High  Expansion  Nozzle 

Ar  indicated  earlier,  the  effect  of  viscosity  on  high  expansion  ratio  nozzles  with 
swirl  are  considerably  greater  than  tnat  on  C-D  nozzles.  Supersonic  flows  through 
a  contoured  nozzle  with  an  expansion  ratio  of  272  :  1  as  that  given  in  Chapter  3 
were  computed  by  using  the  PNS-AD!  algorithm  (Eq.  (3.40)).  A  75  x  50  -rid  as 
shown  in  Fig.  30  was  dsed  and  the  Reynolds  number  based  on  inlet  (throat)  radius 
and  inflow  conditions  was  1.4  x  105  for  both  laminar  and  turbulent  ca'cuiMionii. 
The  convergence  rates  for  zero-swirl  and  typical  swirling  case  are  shown  in  Fig.  65. 
I  lie  results  show  that  swirling  has  minor  effects  on  the  convergence  rate  for  the 
1’NS-ADI  algorithm. 
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The  Mach  number  contours  from  laminar  results  for  5,  =  0  and  S,  =  0.521 
are  shown  In  Fig.  69.  The  presence  of  swirl  increases  the  axial  velocity  near  the 
centerline  and  thus  results  in  a  shifting  of  iso-Mach  lines  and  the  weak  oblique 
shock.  The  resulting  boundary  layers  are  thicker  in  these  calculations  than  those 
in  C-D  nozzles  even  for  a  Reynolds  number  as  high  as  1.4  x  10®.  A  much  thicker 
boundary  layer  can  be  expected  for  lower  Reynolds  number  flows  where  inviscid 
assumptions  appear  to  be  inadequate. 

Turbulent  results  by  using  Baldwin  and  Lomax  model  |12|  are  shown  in  Fig.  70 
based  on  the  same  Reynolds  number  and  two  inlet  swirl  numbers  of  0  and  0.55. 
Comparing  with  Fig.  69,  relatively  thicker  boundary  layer  is  seen.  Figure  71  plots 
Cp  and  tj„,  as  functions  of  5,  for  both  laminar  and  turbulent  results.  Large  reduc¬ 
tion  in  Cp  and  rjv t  for  both  laminar  and  turbulent  results  can  be  observed.  For  a 
highly  swirled  flow  (5,  =  2.5),  Cp  and  tjv ,  are  about  20%  less  than  those  of  the 
zero-swirl  flow,  even  for  a  moderate  swirl,  a  10%  reduction  in  Cp  and  nvt  may 
occur.  Slightly  less  reduction  in  Cp  and  r)vt  for  turbulent  results  are  noted.  These 
results  demonstrate  the  effects  of  swirl  on  high  expansion  ratio  nozzles  are  much 
more  prominent  than  those  on  mild  C-D  nozzles. 

5.3. 4  Verification  of  Global  Conservation 

To  validate  the  numerical  algorithms,  the  mass  flow  rate  at  each  axial  location 
is  calculated  in  the  analysis  codes.  This  provides  a  back-to-back  chock  for  global 
mass  conservation.  For  the  transonic  results  presented  above,  the  maximum  mass 
error  has  been  maintained  below  1%.  For  the  more  difficult  high  expansion  nozzle 
case,  which  has  the  largest  mass  error  to  date,  the  maximum  deviation  is  about 
0.8%.  This  again  verifies  the  necessity  of  fully  conservative  form  for  the  internal 
flow  calculations. 


CHAPTER  6 


THREE-DIMENSIONAL  NOZZLE  FLOWS 


The  hybrid  upwind-central  algorithms  proposed  in  Chapter  2  are  extended 
to  three-dimensional  viscous  supersonic  calculations  in  this  chapter.  The  three- 
dimensional  thin-layer  Navier-Stokes  equations  are  simplified  by  neglecting  the 
streamwise  diffusion  while  retaining  all  viscous  terms  on  the  cross-stream  plane. 
Both  the  Parabolized  Navier-Stokes  procedure  and  the  time-iterative  TLNS  algo¬ 
rithm  are  studied  for  three-dimensional  nozzle  flowfield  predictions.  These  algo¬ 
rithms  are  formulated  based  on  the  DDADI  splitting  for  the  streamwise  derivative 
and  central  differencing  in  cross-plane  derivatives.  Supersonic  flows  through  a  three- 
dimensional  nozzle  with  a  rectangular  cross  section  are  computed  for  demonstration. 

6.1  Governing  Equations 

The  three-dimensional  Navier-Stokes  equations  in  a  Cartesian  coordinate  sys¬ 
tem  can  be  written  in  vector  notation  as 

dQ  dE  dF  dG  dEv  dFv  0GV  ,  , 

"TT  +  7J—  +  -r-  +  =  ~z —  +  ~r —  +  ~ —  (6.1) 

at  ax  dy  dz  dx  dy  dz 

where  the  dependent  variable  Q  is 


Q  =  (p,pu,pv)pw,e)T . 
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The  inviscid  flux  vectors  E ,  F ,  and  G  are 
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Viscous  terms  are  included  in  the  flux  vectors  Ev,  Fvt  and  Gv  as 
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To  facilitate  computations  on  arbitrary  grids,  the  Cartesian  coordinates  x,  y, 
and  r  are  transformed  to  genera!  coordinates  £,  77 ,  and  f  according  to 

£  =  £(*,y.*) 

r)~ri[x,y,z)  (6.2) 

f  =  f(*.y.*)- 

By  using  the  transformation  defined  by  Eq.  (6.2),  Eq.  (6.1)  can  be  transformed  to 
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The  dependent  variable  now  becomes  Q ,  which  is  defined  by 

t 

Q  =  j{p,pu,pv,pw,t)T 

where  J  is  the  Jacobian  of  the  coordinate  transformation  and  can  be  expressed  as 
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in  which  the  contravariant  velocities  U,  V ,  and  W  in  three  dimensions  are 

U  =  (gU  +  £yV  +  £,ti> 

•r  =  VzU  +  VyV  +  HrW 


W  =  fjU  +  fyV  +  fZU>. 

The  transformed  viscous  flux  vectors  are  defined  by 

Ev  —  {ZzEv  +  T 
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As  mentioned  in  Chapter  3,  the  streamwise  diffusion  can  be  neglected  without 
losing  accuracy  even  for  a  fairly  low  Reynolds  number.  For  three-dimensional  Hows 
inside  a  nozzle,  the  viscous  terms  in  rj  and  f  directions  cannot  be  neglected.  This 
results  in  the  TLNS  equations  in  three  dimensions  as 
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The  viscous  terms  on  the  right  hand  side  can  be  further  rearranged  to 
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where  Fn  and  G(  contain  pure  second  order  derivatives  with  respect  to  rj  and  <;, 
respectively.  All  cross  derivatives  are  included  in  Ff  and  Gn.  The  viscous  flux 
vector  Fn  can  be  expressed  as 


Qlf^  +  +  <*3J% 

Ctj  -+■  Ctj 

d  u  ,  .  dv  i  _  9u> 


Fr,  = 


ft„l£x  ql.-q.ji  1x1  x.  2*-,?jh  t>xl  +  aa-fltt  1ml 

p  ^  7  df)  T  7  dW  ^  7  dn 


1  ^  gut>  ,  _  guu;  1  _  8vw 

+a2*5T  +ft3V  +a*~5T 


with 


»o  — 


-r*(f)?  +  ,i2+T?) 

7cv 


Q2  = 


H  *>»  na 
J  3 


<*3=5^’ 


as  = 


7  3 


ai  =  jihl  +  n]  +  *7*) ) 

*«-5(*».a  +  K +  *».’). 

a0  =  5^*  +  ''y  +  5^*)- 

The  vector  Gf  can  be  obtained  by  replacing  all  77’s  in  the  expression  of  F„  with  ?'s. 
The  vectors  Ff  and  C7?  are 
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+  7  5^  +  18 

.  <3  y  .  _  (3  v  .  c3  u» 

at  +  Jn 

vl<  .i.  31-rJu.  3«I  ,  .lijrJu  Jvl  +  1? 

9rj  <>  2  dn  ’  2  On  2  0  1 

+  'T4u|^  +  12V%%  +  1T7U§*  +  T3U>f* 
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with 


*72  =  5  (^yf*  —  3  ^xCy ) »  "73  =  ~  §  *?*?«)* 

*>4  =  5(’l*fv  ~  §^vf*)i  '*78  =  5(n«f»  +  l^vfy  + 

nro  =  {H’J'fy  -  I’M*)’  "77  =  -  f'f.C*). 

'll  =  5(w*  ~  I'My).  *7®  =  +  5»7.f.)- 

With  the  substitution  of  the  new  expressions  of  Viscous  terms,  Eq.  (6.4)  be¬ 


comes 


9Q  |  dE  |  dF  |  dG  dA, 


at 


a^  acn  a^ 

4-  •r —  4-  *r—  =  ~r —  4-  — —  4-  — —  4-  — — .  (6.5) 

a^  di)  af  a^  aq  a?  a^ 


0-2  Three-Dimensional  Supersonic  Algorithms 

Numerical  algorithms  for  the  solution  of  Eq.  (6.5)  can  be  formulated  in  a  num¬ 
ber  of  ways.  Based  upon  the  results  from  Chapter  3,  the  algorithm  for  three- 
dimensional  flows  will  be  formulated  according  to  flux-vector  splitting  in  the  stream- 
wise  direction  and  central-difTerencing  in  cross-stream  directions.  Before  discussing 
the  details  of  numerical  algorithms  for  the  vector  governing  equations,  the  Fourier 
stability  analysis  for  a  scalar  modeled  equation  is  studied. 


6*2.1  Stability  Analysis  of  the  Scalar  Equation 


The  three-dimensional  Burger's  equation, 


du  ,  du 

T,+a  5J+0' 


du  du  du  d2u 

^+t’Ty+ClT^-,‘(JS 


dx 


d2u, 
dz7  ‘ 


(6.6) 


is  chosen  as  the  modeled  equation.  This  modeled  equation  implies  only  the  stream- 
wise  (x)  direction  is  flux-vector  split,  while  the  remaining  derivatives  on  the  cross 
plane  (y  and  z  directions)  are  to  be  evaluated  according  to  central  differences. 

For  simplicity,  we  restrict  only  to  first-order  upwind  differencing  for  the  dis¬ 
cussion  of  stability  analysis  although  the  numerical  computations  shown  later  are 
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baaed  on  second-order  upwind  differencing.  By  using  the  line  Gauss-Seidel  version 

i 

for  the  DDADI  method,  the  discretized  equation  for  Eq.  (6.6)  can  be  expressed  as 
a  forward  marching, 

-  *v> + -  -*'*  m 

and  a  backward  marching, 

\d  +  A t(bj-  -  *5^)  +  -  M^yllAu  =  -  Air'.  (6.8) 


dz 


In  Eqs.  (6.7)  and  (6.8),  the  right  hand  side  residuals  rn  and  r*  are 

r"  =  +  |0-  ^  +  haJt  +  caJt  _  ^  + 

Ax  1  dx  dy  dz  dz'1. 


"  u'*.;.k  ,  f  +5u  .  ,  ^  .  du  ,d2u  d2u,,. 

r  =  a  - - - - —  +  o  —  +6  -+c - u  -  H - ) 

A  z  1  dx  dy  dz  ^dy*  dz')] 


and  the  quantity  d  is  the  diagonal  element  defined  by 

d  =  1  +  — (a+  -  a" ). 

Ai 

Equations  (6.7)  and  (6.8)  are  based  on  a  straightforward  extension  of  the  two- 
dimensional  algorithm.  This  implies  that  the  y  and  z  derivatives  (from  both  inviscid 
and  viscous  terms)  are  treated  implicitly.  Consequently,  the  resulting  left  hand  side 
matrices  of  Eqs.  (6.7)  and  (6.8)  are  very  expansive  to  solve  due  to  their  high  band¬ 
width  structure. 

A  more  practical  way  to  solve  these  two  equations  is  to  factorize  the  left  hand 
side  operators  of  Eqs.  (6.7)  and  (6.8).  This  results  in  the  factored  forward  marching, 

c?  d  ^2 

\d  +  At(6— -  -  ii—yd-'ld  +  Af(c—  -  mt-t)|Au*  -  -Afr"  (6.9) 
dy  dy1  dz  dz 21  v  ’ 
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and  the  factored  backward  marching, 

\d  +  At(i>^  ~  +  At(c^  ~  ^^j)lAu  ®  “A<r*-  (610) 

The  solutions  of  these  factored  equations  require  alternating  sweeps  in  y  and 
2  directions,  each  sweep  involves  only  consecutive  solutions  of  a  scalar  tri-diagonal 
matrix  for  the  present  scalar  modeled  equation.  In  the  vector  governing  equations  or 
interest,  the  factored  algorithm  requires  the  solutions  of  a  block  tri-diagonal  matrix 
with  a  block  size  of  5  x  5. 


The  amplification  factor  for  the  unfactored  forward  marching  Eq.  (6.7)  is  then 


where,  £>*  and  C’  are 


D'u  =  1  -  o~  (cos  ioz  +  i  sin  wr) 
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and  and  u,  are  von  Neumann  numbers  defined  by 


nAt  nAt 


The  wavenumbers  in  x,  y,  and  z  directions  are  represented  by  wt,  u>v,  and  w,, 
respectively.  The  overall  amplification  factor  for  the  unfactored  method  is  then 

Su  —  9 u9u  • 

The  amplification  factors  for  the  factored  forward  and  backward  marching  are 

_  d:  +  c, 

S/  C;  +  c, 

and 

..  =  py  +  Cf 
9/  Cy  +  C/  ’ 

respectively.  The  quantity  Cj  is  due  to  the  approximate  factorization  and  is 

C/  =  (1  +  az  -  °x)V°y  oinu»v  +  2uv(l  -  cos  wv)][t<7,  sinw,  +  2^,(1  —  cos  w,)|. 

Th*  overall  amplification  factor  for  the  factored  method  is  gj  =  9)9)' ■ 

Similar  to  the  discussion  in  Section  3.3.5,  the  stability  results  are  presented  for 
two  special  cases,  they  are 

1.  subsonic  \o+  =  -a~  —  cv  =  az  —  vy  =  vt  —  CFL 

2.  supersonic:*?*  =  ov  =  <?f  =  uy  =  =  CFL .  o~  =  0. 

Here,  the  first  case  simulates  the  subsonic  flow,  while  the  second  case  is  analogous 
to  the  supersonic  flow. 

The  results  for  the  first  case  with  a  CFL  number  of  10  are  shown  in  Fig.  72 
and  Fig.  73  for  the  unfactored  and  factored  schemes,  respectively.  These  two  figures 
plot  the  amplification  factors  versus  the  wavenumber  uy  and  u for  three  typical 
x-direction  wavenumbers,  wz  =  0,  u/z  -  tt/2,  and  uz  =  rr.  As  can  be  seen,  the 
approximate  factorization  of  the  factored  scheme  results  in  a  much  higher  ampli¬ 
fication  factor  near  the  high  wavenumber  region  as  compared  to  the  unfactored 
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case.  Therefore,  the  factored  scheme  is  expected  to  give  much  slower  convergence 
than  the  unfactored  scheme.  The  stability  results  for  the  second  case  are  shown  in 
Fig.  74  and  Fig.  75.  Similar  effects  of  approximate  factorization  as  in  the  first  case 
are  clearly  shown. 

The  results  above  indicate  that  both  the  factor  and  unfactored  three* 
dimensional  DDADI  schemes  are  unconditionally  stable  for  the  three-dimensional 
Burger’s  equation  Furthermore,  the  necessary  approximate  factorization  for  a  prac* 
tical  three-dimensional  DDADI  algorithm  results  in  increasing  the  eigenvalues  near 
the  high  wavenumber  region.  This  implies  that  convergence  of  the  three-dimensional 
DDADI  algorithm  is  inferior  to  that  of  the  two-dimensional  DDADI  algorithm  given 
in  Section  3.4.  In  fact,  due  to  similar  approximate  factorization  required,  we  ex¬ 
pect  this  three-dimensional  algorithm  will  give  similar  convergence  as  that  of  the 
two-dimensional  central-differenced  ADI  scheme. 

0-2.2  Numerical  Algorithms  of  the  TLNS  Equation 

Similar  to  the  formulations  of  the  DDADI  algorithm  in  axisymmetric  two  di¬ 
mensional  flows,  we  first  split  the  streamwise  flux  vector  E  according  to 

E  =  E+  +  £-. 

This  splitting  can  be  done  by  using  either  the  characteristic  splitting, 

E+  =  [M(\t  M"l)Q  E~  =  (M<A'A  r~')Q, 
or  the  pressure  gradient  splitting, 


pU 

r  0 

puU  -f  u£zp 

■  *'  =  7 

(1  -  u>)(xp 

PVU  + 

(1  -  U>)tyP 

pvuU  +  u£tp 

(1  -  u){rp 

.  (e  +  P)U  - 
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Figure  75.  Stability  results  of  3-D  Burger  s  equalii 
for  supersonic  case 
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The  matrices  A/{  and  M ^  1  are  the  right  and  the  left  eigenmatrices  of  the  Jacobian 
matrix  A  ~  0E/3Q ,  and  the  diagonal  matrices  A,  and  A,  are 

\  v 

+  _  +  |Af|  A(  -  |A*| 

A( - - - ,  A{  =  - - - 

where  A,  =  iiag[U,U,U,U  +  C, ,  V  -  C()  with  Cf  =  Jt\  +  {J  +  (|e.  The  pa- 
rarneter  u>  of  the  pressure  gradient  splitting  in  three  dimensions  is  equal  to  unity  if 
U  >C(  and  must  satisfy 

w  <  1  +  h  -  1  )M\ 

if  U  <  C(,  where  Af(  is  the  streamwise  Mach  number  (A/(  =  U/'C {). 

The  flux  vectors  £“*,  F,  and  G  in  Eq.  (6.5)  can  be  linearized  by  the  truncated 
local  Taylor  series, 

(£±)nM  =  [E±)n  +  A±AQ 


Fn  +  l  =  Fn  +  BAO 
Gn+l  =  Gn  +  CAQ 


in  which,  A  ,  B,  and  C  are  Jacobian  matrices  of  £'~,  1' ,  and  G,  respectively. 
Viscous  terms  containing  pure  second  order  derivatives  can  be  linearized  according 
to 

where  the  viscous  Jacobians  Bv  and  Cv  can  be  expressed  in  a  similar  way  as  in 
Eq.  (3.10)  For  example,  Bv  is 


Sr  =  R^ 


d  ,  dQ 


dr ;  dQ 


d  dQ2 , 


dn  dQ  ' 
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where  the  matrices  R\  anti  i?2  are 


r  0 

0 

0  o  o  ■ 

0 

Oi 

Q  2  Q3  0 

Ri  = 

0 

q2 

a<  as  0 

0 

03 

Q5  Oc  0 

.03 

0 

0 

D  a5. 

■  0 

0 

0 

0 

0  ■ 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

.  Oo 

a  i  —  a 

u 

a*  -  au 

”Q(i 

<»2  - 

3 

2 

2 

and  the  vectors  (?|  and  Q?  are 

Qi  =  (uv,u,v,  ,w,vw)T,  Q2  =  (e/p,u7,v^,w7,uv)T. 


The  matrix  C v  can  be  obtained  by  replacing  all  a’s  with  0's  in  the  expressions  of 
Ri  and  R 2  above.  The  cross-derivative  viscous  flux  vectors  in  Eq.  (6.5)  will  not  be 
linearized  because  the  linearization  of  these  vectors  will  results  in  a  high  band-width 
left  hand  side  matrix. 

Direct  application  of  the  DDADI  splitting  to  the  £  derivative  of  Eq.  (6.5)  results 
in  the  unfactored  forward  marching 


I D’  +  -r  —  C  -  —  Bv  -  —  CV)|AQ*  =  -At/?n 

dr]  d <;  d<; 


(6.11) 


and  the  unfactored  backward  marching 


()  d  d  d 

1 D'  r-  +  —  C  -  B v  -  —  CV)|AQ  =  -Ar/r 

di]  di -  di]  0$ 


(6.12) 


where  the  residuals  Rr'  and  R’  are 


n  (^r-2 (jtxj.kV  +  (E+.2.,.kv 


Rn  = 


As 


2A£ 


-l-  ■ 


,d£- 


01  Of,  0Fn  OF -  OGrj  d  G’; 


d£  dr ;  0>]  Or]  d f  d( ‘ 
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and 


P.  (£>.., .*)"+'  -  .(W  -  +  <*,«.,.»> 

* - A? 


n  +  l 


-  /c- 


i££l  ££  ac  _  af, 

+  +  dr?  +  df  dr)  drj 

The  diagonal  matrix  D'  is 


2A£ 

acn  _  dG< 
d<;  dc  1  ’ 


D'  =  I  +  (l  +  k/2)^[A+  -  A~) 

where  the  quantity  ic  is  0  for  first  order  and  is  i  for  second  order  upwind  differencing 
in  C 

To  avoid  the  solution  of  a  high  band-width  matrix,  the  left  hand  side  operators 
of  Eq.  (6.11)  and  Eq.  (6.12)  must  be  further  factorized  into  two  one-dimensional 
operators.  This  results  in  the  factored  forward  marching, 

i D'  +  Al ( Jjfl  -  -BMDT'ID'  +  Al(  Jtc  -  ^C„)|AQ'  =  -Al ft"  (C.13) 

and  the  factored  backward  marching, 


( D'  +  A t(4~B  -  ~Bv)}(DYl\D'  +  A t(~C  -  |-C,)|AQ  =  -Af/T.  (6.M) 

OT)  dr)  d<;  d<; 


Note  that  the  corresponding  Parabolized  Navier-Stokes  procedure  can  be  ob¬ 
tained  by  neglecting  the  flux  vector  E~  and  its  Jacobian  A~  in  Eq.  (6.13).  Either 
the  traditional  presaure  gradient  splitting  or  the  characteristic-based  splitting  ran 
be  chosen  to  form  thio  three-dimensional  PNS  algorithm.  Also,  'lie  combination  of 
Eq.  (G.13)  and  Eq.  (6.14)  provides  a  three-dimensional  TENS  solver. 

To  assess  the  numerical  efficiency  of  these  three-dimensional  algorithms,  the 
supersonic  flow  through  a  15°  expanding  three-dimensional  nozzle  with  rectangular 
cross-sections  was  chosen  for  numericai  experiments.  The  nozzle  geometry  is  shown 
in  Fig.  76  and  the  grids  on  the  inlet  plane  and  the  side  wall  surface  are  shown  in 
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Fig.  77.  A  constant  Mach  number  of  1.2  and  zero  contravariant  ve!ot..ies  of  V  and  IV 
were  imposed  at  the  inlet.  The  flow  was  assumed  laminar  and  the  Reynolds  number 
was  taken  to  be  10s  based  on  the  inlet  hydraulic  radius  and  inflow  properties. 

The  numerical  experiments  were  done  for  both  the  three-dimensional  PNS  and 
TLNS  algorithms  with  approximate  factorization.  Typical  convergence  curves  for 
the  PNS  algorithm  are  shown  in  Fig.  78.  As  we  can  see,  due  to  the  additional 
approximate  factorization  of  the  left  hand  side  operator,  the  quadratic  convergence 
in  two-dimensional  PNS  procedures  cannot  be  obtained  for  the  three-dimensional 
PNS  solver.  The  optimum  CFL  number  for  this  .  sc  is  20  and  this  results  in  300 
local  iterations  to  reach  7  orders  of  magnitude  reduction  in  the  L-2  norm.  But 
acceptable  convergence  (5  orders  reduction  in  the  L-2  norm)  can  be  achieved  in  40 
iterations.  However,  this  time-iterative  three-dimensional  PNS  algorithm  has  been 
found  to  be  very  robust  and  is  insensitive  to  grid-stretching  in  the  {  direction. 

The  convergence  for  the  three-dimensional  TLNS  procedure  is  shown  in  F:g.  79 
for  a  optimum  CFL  number  of  20.  The  initial  conditions  for  this  calculation  were 
obtained  from  the  corresponding  converged  (5  orders  of  magnitude  reduction  in 
the  L-2  norm  at  each  cross-plane)  PNS  solutions.  As  is  seen,  for  300  iterations, 
the  L-2  norm  drops  7  orders  of  magnitude,  which  is  about  the  same  rate  as  a 
two-dimensional  central-differenced  ADI  solver,  as  we  predicted  from  the  stability 
analysis  given  in  Section  6.2.1. 

0.3  Flowfield  Predictions 

The  test  case  for  three-dimensional  flowfield  predictions  wa;i  the  supersonic 
flow  through  a  three-dimensional  nozzle  with  rectangular  cross-sections  is  shown  in 
Fig.  80.  The  wall  contour  of  this  nozzle  was  chosen  to  be  the  same  as  that  of  the 
272  :  1  axisymmetric  contoured  nozzle  previously  given.  This  nozzle  has  a  constant 
width  of  30  mm.  in  the  y  direction.  Due  <.o  the  symmetry  conditions,  only  one 


A-A  section 


B-B  section 


Figure  77.  5  x  30  x  30  grid  of  15°  expanding  3-D  nozzle  for  conv-rgence  test 


rnetry  of  3-D  flowfield  prediction 
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quadrant  on  the  cross-section  need  to  be  calculated. 

t 

The  75  x  30  x  30  grid  system  with  75  in  the  x-direction  and  30  x  30  on  one 
quadrant  of  the  cross-plane  is  shown  in  Figs.  81-82.  Figure  81  shows  the  grid  on 
the  side  wall  of  the  nozzle.  As  is  seen,  a  strong  clustering  near  the  top  is  chosen  to 
resolve  the  boundary  layer.  Typical  grids  on  the  cross-planes  are  shown  in  Fig.  82 
for  both  the  inlet  and  exit  planes.  This  figure  also  shows  strong  stretching  near  the 
side  wall  and  the  top  due  to  the  boundary  layers  in  the  y  and  z  directions. 

The  inlet  Mach  number  was  assumed  to  be  uniformly  1.02  and  the  gas  prop¬ 
erties  described  in  Section  3.6  were  imposed  at  the  inlet.  This  resulted  in  a  nozzle 
Reynolds  number  based  on  the  throat  hydraulic  radius  of  1.5  x  101  and  laminar 
flow  was  assumed.  The  PNS  procedure,  Eq.  (6.13),  was  utilized  to  perform  the 
calculation. 

The  flowfield  results  are  shown  at  several  locations  indicated  in  Fig.  83.  The 
computed  Mach  number  contours  at  locations  A  and  B  are  shown  in  Fig.  84.  This 
figure  shows  quite  different  results  from  those  of  the  axisymmetric  calculations  pre¬ 
sented  in  previous  chapters  due  to  the  three-dimensionality.  Although  not  shown 
here,  these  Mach  number  contours  are  more  similar  to  corresponding  planar  two- 
dimensional  results.  The  wiggles  of  the  contours  near  the  exit  and  the  center  plan*- 
are  possibly  due  to  insufficient  grid  resolution  and  the  reflection  of  a  weak  oblique 
shock  at  the  center  plane. 

The  slreamwise  velocity  contours  on  several  cross-planes  (indicated  in  Fig.  83 
as  locations  C,  D,  E ,  and  F)  are  shown  in  Figs.  85-87.  In  these  figures,  the  growing 
of  the  boundary  layer  thickness  near  the  side  wall  and  the  top  is  clearly  shown 
Figure-  88-90  show  the  cross-stream  velocity  vector  plots  at  locations  C,  D,  E ,  and 
F .  Secondary  flow  patterns  and  the  development  of  vortices  of  the  three-dimensional 
boundary  layer  near  the  side  wall  and  the  top  are  observed.  The  secondary  (low 
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Figure  81.  75  x  30  surface  grid  on  the  side  wall  of  the  3-D  nozzle 
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C-C  section  exit 


B-B  section  inlet 


Figure  82.  30  x  30  grids  on  the  inlet  and  exit  cross-plane. 


Figure  84.  Mach  number  contours  at  locations  A  and  B 
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patterns  adjacent  to  the  side  walls  a a  shown  in  Figs,  ft 9  and  00  also  explains  why 
the  boundary  layer  along  the  side  wall  is  thicker  near  the  center  plane. 


CHAPTER  7 


SUMMARY 


Implicit  time-dependent  scheme*  have  been  successfully  applied  to  solve  the 
compressible  thin-layer  Navier-Stokes  equations  in  multi-dimensions.  Preliminary 
applications  of  the  implicit  algorithm  to  the  one-dimensional  Euler  equations  were 
studied  by  using  spatial  discretizations  based  on  both  central-difTerencing  and  flux- 
vector  splitting  upwind-differencing.  Both  differencing  methods  were  shown  to  give 
rapid  convergence  and  accurate  solutions.  The  Fourier  stability  analysis  has  been 
studied  for  either  differencing  method.  The  results  were  shown  to  provide  useful 
information  about  the  convergence  criteria.  In  particular,  the  explicit-like  CFL  lim¬ 
itation  of  the  one-dimensional  upwind  scheme  when  using  approximate  Jacobians 
was  successfully  predicted  from  stability  analysis  and  later  on  confirmed  by  nu¬ 
merical  experiments.  The  preparatory  investigations  on  one-dirnensional  flows  also 
provide  informative  results  that  are  extendible  to  multi-dimensions. 

For  two-dimensional  calculations,  the  ADI  scheme  based  on  central-differencing 
was  formulated  to  solve  the  TLNS  equations  in  a  cylindrical  coordinate  system. 
The  effectiveness  of  this  ADI  scheme  was  tested  by  calculating  typical  subsonic, 
transonic,  and  supersonic  flows  through  nozzles.  The  results  showed  that  the  con¬ 
vergence  rates  of  subsonic  and  transonic  cases  are  slow  (when  compared  to  one- 
dimcnsional  calculations)  but  acceptable.  When  the  flow  is  predominantly  super¬ 
sonic,  the  ADI  scheme  has  proven  to  be  inefficient  and  sometimes  even  unstable. 

Based  upon  the  physical  character  of  viscous  supersonic  flows,  a  hybrid  dis¬ 
cretization  composed  of  central  differencing  in  the  streamwise  direction  and  second- 
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order  upwinding  in  the  cross-stream  direction  wu  proposed.  Stability  analyses  on 
a  modeled' equation  were  considered  for  the  fully  implicit  and  three  approximate 
factorisation  procedures  based  on  this  hybrid  discretization  scheme.  The  results 
showed  that  all  four  algorithms  are  unconditionally  stable  for  the  Burger's  equation- 
Further,  the  line-relaxation  version  of  the  DDADI  algorithm  gives  the  eigenvalues 
of  the  amplification  matrix  that  approach  the  fully  implicit  limit. 

On  the  basis  of  encouraging  stability  results,  four  algorithms  indicated  above 
were  then  applied  to  solve  the  TLNS  equations  for  flows  through  nozzles.  Of  the 
three  approximate  techniques,  the  DDADI  scheme  suggested  by  Lombard  [19|  is 
shown  to  require  the  least  number  of  iterations,  but  in  terms  of  CPU  time,  the 
PNS-ADI  scheme  developed  in  this  study  is  as  fast  as  the  DDADI  scheme.  The 
standard  ADI  factorization  arising  from  traditional  ADI  schemes  18,9,10)  proves  to 
be  the  most  inefficient  in  terms  of  both  the  number  of  iterations  and  the  CPU  time 
required.  The  direct  method  devised  from  the  physics  of  high  Reynolds  number 
viscous  flows  seems  to  be  particularly  suited  for  supersonic  problems,  but  rapid 
convergence  of  the  DDADI  and  the  PNS-ADI  schemes  allows  them  to  surpass  the 
direct  method  in  terms  of  CPU  time  required.  However,  the  direct  method  has 
proven  to  be  more  robust  than  any  approximate  schemes  by  noting  that  the  CFL 
number  can  be  as  high  as  10l°  without  losing  stability.  Numerical  experiments 
regarding  boundary  conditions  and  Jacobian  matrices  have  indicated  that  implicit 
boundary  conditions  toge  er  with  true  Jacobians  of  the  split  flux  play  a  decisive 
role  on  convergence. 

The  solutions  of  present  upwind-central  differencing  algorithms  were  compared 
to  those  of  MOC  calculations,  excellent  agreements  on  the  wall  pressure  distribution 
and  Mach  number  contours  were  demonstrated. 
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For  the  first  time,  proper  downstream  boundary  conditions  are  applied  for 
the  subsonic  portion  of  the  outflow.  These  downstream  boundary  conditions  were 
shown  to  be  capable  of  allowing  supersonic  solutions  to  respond  to  nozzle  back 
pressure  conditions  as  they  should  do  in  realistic  situation.  The  calculations  with 
recirculation  and  reentry  flows  at  the  exit  plane  caused  no  difficulty,  and  the  results 
showed  that  different  exit  pressures  would  alter  the  nozzle  boundary  layer  charac¬ 
teristics  near  the  exit  plane.  The  extrapolation  conditions  that  are  normally  used 
were  shown  to  give  solutions  corresponding  to  one  specific  back  pressure  condition. 
A  series  of  results  showing  the  effects  of  variations  in  back  pressures,  wall  tempera¬ 
tures,  and  nozzle  Reynolds  numbers  are  given  for  both  a  conical  nozzle  and  a  high 
expansion  ratio  contoured  nozzle.  The  effects  of  turbulence  on  supersonic  nozzle 
flows  with  separation  were  investigated  by  solving  the  Reynolds  averaged  Navier- 
Stokes  equations  with  the  Baldwin  and  Lomax  model.  The  global  characteristics  of 
turbulent  flows  were  properly  resolved  by  using  this  algebraic  turbulence  model. 

The  results  of  testing  on  the  global  mass  conservation  indicate  that  global 
mass  errors  can  be  kept  below  1%  when  fully  conservative  form  is  used,  while  quasi 
conservative  form  may  give  a  global  error  as  large  as  30%  even  in  flowfields  without 
discontinuities. 

Along  with  the  development  of  Navier-Stokes  algorithms,  the  parabolized  pro¬ 
cedures  were  also  investigated.  Unlike  the  traditional  approach,  the  PNS  procedure 
was  obtained  from  the  Lime-dependent  general  flux-vector  split  TLNS  equations, 
for  which  the  streams  'se  flux  vector  has  been  split  into  two  parts  corresponding 
to  downstream  and  upstream  characteristics.  By  omitting  the  parts  with  upstream 
characteristics,  the  whole  equation  set  was  made  parabolic  in  the  slreamwise  di¬ 
rection.  With  this  approach,  a  distinct  PNS  formulation  can  be  obtained  for  each 
type  of  flux-vector  splitting  considered,  Two  examples  were  chosen  for  demonstra- 
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tion.  The  traditional  PNS  formulation  is  obtained  by  using  a  pressure  gradient 
splitting.  iThe  use  of  characteristic-baaed  flux  splitting  yields  a  PNS  algorithm  that 
includes  only  the  downstream  characteristics.  Stability  results  showed  that  this 
characteristic-based  PNS  algorithm  is  stable  for  space-marching  and  numerical  re¬ 
sults  indicated  that  it  provides  solutions  that  are  identical  to  the  classical  pressure 
gradient  split  PNS  formulation  and  in  excellent  agreement  with  the  TLNS  solutions. 

One  advantage  of  the  present  PNS  algorithm  when  compared  to  non-iterative 
space-marching  procedures  is  that  the  current  approach  requires  no  safety  factor. 
Comparisons  of  the  pressure  gradient  splitting  PNS  calculations  with  the  TLNS 
solutions  show  that  the  introduction  of  a  safety  factor  deteriorates  the  solution 
accuracy.  The  necessary  local  iterations  on  each  £  plane  for  the  time-iterative  PNS 
algorithm  result  in  more  computational  time  than  the  traditional  non-iterative  PNS 
procedure.  However,  it  has  been  shown  that  this  CPU  time  deficit  is  partially  offset 
in  that  the  local  iterations  allow  the  £  derivative  to  be  formulated  in  a  conservative 
form  so  that  variable  step  sizes  in  £  can  be  used. 

The  global  pressure  iteration  procedure  in  the  traditional  PNS  approaches  has 
been  interpreted  as  a  TLNS  procedure.  These  mathematically  and  physically  well- 
posed  TLNS  procedures  based  on  approximate  factorization  are  sugjp.'aisd  instead  of 
the  traditional  global  pressure  iteration  procedures  based  on  an  arbitrary  relaxation 
of  the  pressure  field. 

Numerical  algorithms  for  computing  viscous  swirling  nozzle  flows  have  also  l>crn 
studied  by  using  time-iterative  implicit  schemes.  The  implicit  ADI  and  the  PNS- 
ADI  algorithms  are  utilized  to  solve  transonic  and  supersonic  swirling  flows,  respec¬ 
tively.  These  algorithms  prove  to  be  equally  efficient  for  swirling  two-dimensional 
calculations.  The  combined  effects  of  viscosity  and  swirling  on  the  flowfield  and 
the  integral  nozzle  performance  are  investigated  for  transonic  and  supersonic  Hows 
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through  mild  to  high  expansion  ratio  nozzles.  Viscous  calculations  are  performed  for 
three  nozzle  geometries  previously  investigated  by  Dutton  (34|.  These  results  vali¬ 
date  the  inviscid  assumptions  for  high  Reynolds  number  flows  and  show  how  rapidly 
nozzle  performance  deteriorates  with  the  Reynolds  number.  For  the  high  expansion 
ratio  contoured  nozzle  and  the  plug  nozzle,  the  combined  effects  of  swirling  and 
viscosity  have  significant  influence  on  the  flowfielda  and  the  nozzle  performance. 

Finally,  the  algorithms  developed  for  axisymmetric  two-dimensional  flows  are 

extended  to  solve  the  three-dimensional  TLNS  equations.  These  three-dimensional 
* 

algorithms  are  based  upon  the  DDADI  splitting  for  the  streamwisc  flux  vector  and 
an  additional  approximate  factorization  of  the  left  hand  side  operator.  The  optimum 
CFL  number  reduces  to  the  order  of  10  and  it  gives  slow  convergence  due  to  this 
approximate  factorization.  Both  the  PNS  and  TLNS  procedures  are  formulated  for 
three-dimensional  calculations.  Typically,  acceptable  convergence  can  be  achieved 
by  40  local  iterations  for  the  PNS  algorithm.  Although  this  convergence  rate  is  not 
competitive  with  that  of  two-dimensional  algorithms,  these  three-dimensional  algo¬ 
rithms  Drove  to  be  robust  and  are  insensitive  to  grid-stretching  along  the  streamwise 
direction. 

As  a  final  comment,  although  all  the  results  shown  in  this  study  are  compula¬ 
tions  of  internal  flowfields  of  nozzles,  all  numerical  algorithms  developed  here  are 
applicable  to  external  flows  as  well. 
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